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PREFACE 


As far as the laws of mathematics refer to reality, they are not certain; and as 
far as they are certain, they do not refer to reality. - A. Einstein 


Plasmas, such as those found in the space environment or in fusion energy devices, are often 
modeled as electrically conducting fluids. Fluids, and thus plasmas, when energetically 
stirred, have regions of highly nonlinear, chaotic behavior known as turbulence. Although 
turbulence in fluids and plasmas is often the critical factor in determining overall long-term 
behavior, a general theory of turbulence has been elusive. Establishing such a general theory 
is a ‘grand challenge’ to modern science. 

The present work describes a statistical theory concerning a certain class of nonlinear, 
finite dimensional, dynamical models of turbulence. These models arise when the partial 
differential equations describing incompressible, ideal ( i.e ., non-dissipative) homogeneous 
fluid and magnetofluid (i.e., plasma) turbulence are Fourier transformed into a very large 
set of ordinary differential equations. These equations define a divergenceless flow in a high- 
dimensional phase space, which allows for the existence of a Liouville theorem, guaranteeing 
a distribution function based on constants of the motion (integral invariants). The novelty 
of these particular dynamical systems is that there are integral invariants other than the 
energy, and that some of these invariants behave like pseudoscalars under two of the discrete 
symmetry transformations of physics, parity and charge conjugation. 

In this work, we show that the ‘rugged invariants’ of ideal homogeneous turbulence 
are, in fact, the only significant scalar and pseudoscalar invariants. Furthermore, we show 
that the existence of pseudoscalar invariants causes the symmetry of the original equations 
to be dynamically broken, inducing a nonergodic structure on the associated phase space. 
Although realistic turbulence is non-ideal, it also has no existing viable statistical or kinetic 
theory to describe it. In lieu of this, and perhaps as a remote precursor, what is presented 
here is the complement: a statistical mechanics of ideal homogeneous turbulence. This, in 
turn, may bring us closer to meeting the ‘grand challenge.’ 

This work began two decades ago, when the author noticed that canonical ensemble 
predictions did not match the results of numerical experiments by a small but significant 
amount. Following this lead, an interesting story unfolded, leading to an understanding of 
what caused the initial mismatch, and the development of the statistical mechanics of an 
intriguing class of conservative, nonlinear dynamical systems. What is presented here is not 
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an extension of chaos theory, which deals with motion in low-dimensional phase spaces, but 
rather a statistical theory of conservative dynamical systems with relatively large, but finite, 
dimensional phase spaces. 

A note on references: The references given here are mainly the ones that have been useful 
to the author in his research. No attempt has been made to be thorough in the citing of 
references that might be concerned with the general topic of ideal homogeneous turbulence. 
History shows us that many people work in parallel on any given research topic, and that 
essentially identical results are often independently produced. However, the modern age 
has the advantage of electronic search and communication systems, and in using these, the 
author has found no similar work addressing the central subjects presented herein: Broken 
symmetry and nonergodicty in ideal homogeneous turbulence. 

I would like to thank David Montgomery for introducing me to this subject. I thank 
Robert Rubinstein for his review of this document and for catching some errors. Of course, 
any errors that remain are mine alone. (I would appreciate being informed of any errors, 
either typographical or conceptual.) Finally, I would like to thank all the other individuals 
and institutions that have helped support this work over the past two decades. 


John V. Shebalin 

j .shebalin@jsc.nasa.gov 

Advanced Space Propulsion Laboratory 
National Aeronautics & Space Administration 
Lyndon B. Johnson Space Center 
Houston, Texas 77058 
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Chapter 1 

The Problem of Turbulence 


1.1 Introduction 

Turbulence in fluids appears to be a ubiquitous phenomenon, occurring in astrophysical, 
geophysical, and engineering flows. The historical course of fluid analysis, however, was 
to initially concentrate on the examination of laminar flow problems, as exact solutions 
were often possible. This allowed fluid dynamicists and mathematicians to gain analytical 
practice by using tools established in other fields, such as potential and complex variable 
theory, and to build a base from which to gain further understanding of the more complicated 
motions of fluids. Although this practice was useful and necessary, observation of real fluids 
indicates that there are often regions where the flow seems convoluted, unsteady, and prone 
to apparent randomness, involving the generation of fluctuating structures termed ‘eddies’ 
whose instantaneous sizes vary over a wide range. 

As the energy being input into a fluid increases, there is typically a transition from 
laminar to turbulent flow, and this occurrence can be analyzed by assuming that a small 
perturbation to laminar flow occurs, linearizing the relevant equations in terms of this small 
perturbation, and studying its subsequent growth. This step produced and still produces 
valuable insights into fluid and magnetofluid mechanics and establishes a bridge between 
laminar and turbulent flow. In the past, it also allowed research to proceed, as the study of 
fluids could still be approached analytically, with quill and parchment, as it were, in the era 
before the widespread availability of high-speed electronic computing machines. 

Turbulence, the final step and most general state of an energetic fluid, has proved, how- 
ever, to be far more resilient to detailed mathematical analysis than the preceding two steps. 
Even with the advent of modern computing, a full and detailed study is still difficult, if not 
currently impossible, due to the extremely large number of degrees of freedom in any realistic 
turbulent flow, in comparison to the capacity of either available or envisioned computers to 
simulate fluid motion. Computers with their limited memories can only capture an exceed- 
ingly small fraction of the degrees of freedom inherent in a fully turbulent, ‘high Reynolds 
number’ flow. To move beyond this impasse, various analytical theories of turbulence have 
been put forward, based on intelligent guesswork. The most successful of these is the rela- 
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tively simple dimensional analysis of Kolmogorov, leading to the famed energy spectral index 
of -5/3 for the so-called inertial range. Other, more detailed, analytical theories have been 
built on seemingly logical hypotheses, coupled with perturbation techniques, but have not 
been generally successful. 

Turbulence in the physical universe is found in the flow of ionized gases ( i.e ., plasmas 
viewed as magnetofluids) within the sun, the solar wind, and planetary magnetospheres, 
as well as in the fluids found in planetary atmospheres and oceans, and in the myriad 
engineering systems devised by human ingenuity. Turbulence is a flow of energy, from its 
injection at relatively large scales, to its transfer through an inertial range, to its dissipation 
at relatively small scales. This may be thought of as a transfer of energy from larger eddies 
to smaller eddies to smaller ones still, until, at the end of the cascade, flow energy becomes 
simply molecular heat. The challenge is to discover a ‘law of turbulence’ that describes the 
details of this energy transfer. This challenge continues to await the intrepid explorer, who 
may need to invent new theories and mathematical techniques, perhaps building on what is 
known, perhaps making a radical departure. 

The main characteristic of turbulence is that it is a non-equilibrium dynamical process 
in which energy is ultimately dissipated. If energy input is abated at some point, the phe- 
nomenon becomes ‘decaying turbulence,’ and if energy is continually injected into the flow, 
it may be called ‘driven turbulence.’ Realistically, turbulence contains a dissipation mech- 
anism, which, for fluids, is provided by viscosity, and for plasmas, by electrical resistivity 
as well as viscosity. Once heat is produced, thermal conductivity also comes into play, and 
will be important if compressibility is a factor in the dynamics of the fluid. However, since 
density variation is often not a critical factor in turbulent flow, particularly for those of low 
Mach number, then assuming at the outset that turbulent flows are incompressible leads to a 
simpler set of basic equations, while still retaining the essential nonlinear interactions which 
cause turbulence. This approximation will be adopted here. 

Although turbulence involves a large number of interacting degrees of freedom, which 
suggests that statistical mechanics be applied, the presence of dissipation, and the transfer 
of energy generally towards the smaller length scales, requires instead that something like 
kinetic theory be used. Thus, many analytical theories of turbulence have at their heart 
a hierarchy of integral equations, which have been treated by expansions similar to those 
found in quantum field theory, so that such approaches are often called ‘field- theoretic.’ 
In contrast, classical statistical mechanics has conservative systems in equilibrium as its 
domain of applicability, or at least systems that can be approximated in such a way. In the 
real world, it is not possible to turn dissipation off in an astrophysical or geophysical fluid 
(although this does occur in one component of a superfluid). In the world of analytical and 
computational models, however, it is possible to set viscosity and resistivity identically equal 
to zero. Mathematically, this is an example of a singular perturbation problem - setting a 
small parameter multiplying the highest derivative in an equation equal to zero changes the 
fundamental nature of the equation and its solutions. Although this may have its uses, it 
must be remembered that this introduces a layer of approximation that formally disconnects 
real and ideal turbulence. 
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1.2 Ideal Flow 

Neglecting dissipation creates the mathematical entities called ideal fluids and ideal mag- 
netofluids. The advantage gained is that standard equilibrium statistical mechanics can be 
applied (as will be shown), albeit only to model systems. Furthermore, if the flow volume of 
interest is far removed from bounding surfaces, then statistical characteristics can be assumed 
independent of spatial position. Such (real or ideal) turbulence is termed homogeneous, and 
is often also isotropic, unless there is some mechanism imposing a preferred direction on 
the system, such as a constant magnetic field in a magnetofluid, in which case it may be 
anisotropic. Ideal turbulent systems are purely mathematical entities whose governing equa- 
tions contain nonlinear terms identical to those found in the equations of real turbulence. 
Ideal equations lack only the linear dissipative term found in real equations, and though 
this may seem a negligible difference, it is not - it is an essential difference. Nevertheless, 
ideal turbulence has its own attraction and the search for ideal statistical solutions has an 
instructive purpose. 

Since either real or ideal turbulence is highly nonlinear, exact analytical solutions do 
not exist and numerical techniques must be utilized to integrate the equations of motion. 
This introduces further approximation, as computers are finite machines, so that round- 
off errors and finite time-integration steps provide another layer between physical reality 
and computational results. Here, we will use finite Fourier expansions to represent spatial 
variation, which has the benefit of allowing for the exact evaluation of spatial derivatives, 
albeit only down to some minimum wavelength. Although computer models can approximate 
turbulent systems, it must be remembered that their efficacy in representing either real or 
ideal turbulent flows is determined only to the extent that they match physical observations 
in the case of real flows, or some independent statistical theory in the case of ideal flows. 

The purpose of the present work is to develop the statistical mechanics of ideal homoge- 
neous, incompressible, fluid and magnetofluid turbulence. It will turn out that this statistical 
mechanics is an interesting and non-trivial extension of standard canonical ensemble theory, 
due to the presence of significant invariant integrals of the motion, in addition to the en- 
ergy. Some of these additional invariant integrals are pseudoscalars under various symmetry 
transformations of the equations of motion, and this will lead to what has been called broken 
ergodicity. The systems to which this statistical theory is applied are computer models of 
ideal turbulence, and these computer models allow us to run numerical experiments. Since 
these models are highly nonlinear, we do not know in advance, nor can we predict, any pre- 
cise details of the flow, which is highly stochastic. Instead, we can make predictions of ideal 
statistical behavior and compare these with time averages determined through numerical 
integration, thereby testing the statistical theory. 

The development of a statistical theory for ideal turbulence is expedited by the fact that 
the various integral invariants can be expressed as expectation values of quadratic forms, 
which leads to Gaussian distribution functions in the phase spaces under consideration, thus 
allowing a relatively straightforward evaluation of expectation values. Since all fundamental 
interactions in nature are describable by nonlinear partial differential equations, the theory 
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developed here has obvious analogies in other areas of physics and applied mathematics 
where dissipation is either not present or can be ignored. (However, in these other areas, the 
integral invariants are related to forms that may not be quadratic, leading to non-Gaussian 
distribution functions that make expectation values more difficult to evaluate.) Finally, 
there is a remote theoretical possibility that ideal results might be extended into a non- 
ideal domain, and this will be briefly discussed, following the primary item of interest: the 
development of a statistical theory for ideal turbulence. 

1.3 References for Further Reading 

A classic text on turbulence is [Batchelor 53]. 

A more modern text that gives a good overview and critique of the current status of turbu- 
lence research is [Frisch 95]. 

Analytical theories, although they do not provide a solution to the ‘problem of turbulence,’ 
indicate the inherent difficulty of this problem by showing the forests of complication into 
which researchers have been forced to enter. A detailed discussion of these theories, following 
a good preliminary overview of the basics of turbulence, is found in [McComb 95]. 


Chapter 2 

The Equations of Motion 


First, we will establish the basic equations of fluids and magnetofluids, and then reduce 
them to those needed for the study of incompressible ideal homogeneous turbulence. We 
begin by assuming that a continuous, electrically conducting fluid exists, although it is pos- 
sible to start from an underlying set of discrete particles, and use kinetic theory to derive a 
single fluid continuum approximation. In the physics of plasmas, this is a common proce- 
dure, and moves the perspective from microscopic kinetics to macroscopic continua. In the 
continuum approximation, characteristic length scales are much greater than inter-particle 
distances, and charge separation and associated quasi-static electric fields are assumed to 
be dynamically unimportant. However, electrical conductivity is generally still present, so 
that magnetic fields may arise due to self-induced or externally imposed electrical currents. 
The study of such fluids is often termed magnetohydrodynamics (MHD) or magnetofluid 
mechanics , and it contains, as a subset, (nonconducting) fluid mechanics. 

The state of a magnetofluid in a given volume V and interval of time T is completely 
specified if its density p, velocity u, magnetic field B, and pressure p are known at all points 
x of the volume V at all times t within the interval T. Thus, density, velocity, magnetic 
field, and pressure are generally functions of x and t, and may be denoted by p(x, t), u(x, t), 
B(x,t), and p(x,t), respectively. (However, in what follows, one or both of the arguments x 
and t may be omitted for brevity.) There are, of course, other assumptions that can be made, 
resulting in various levels of complexity. For our purpose, the stated assumptions, along 
with incompressibility, will suffice, and the resulting system of equations will be sufficiently 
challenging. 


2.1 The Continuity Equation 


Consider now a small volume 5V within the fluid. If the bounding surface of 5V moves with 
the fluid, then the mass 5m = p5V within 5V is constant: 


dSV 
p a 


+ 6V 


dp 

dt 


d8m 

dt 


5 


0 . 


(2.1) 
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Now, assume that 5V is an infinitesimal cube, which can be represented in Cartesian coor- 
dinates as 5V = 8x5y5z, where 8x — x 2 — x\, etc. Neglecting terms of higher order than 8x, 
we have, using x = dx/dt = u x , 


d8x 

dt 


X2 — X\ 



( 2 . 2 ) 


Since analogous results hold for 8y and 6z, we have 


dSV dSx d8y d8z 

, = oydz , + oxoz , + dxoy , 

dt y dt dt dt 


= 8x8y8z 

\ 

= <5VV • u. 

Combining (2.1) with (2.3) gives 

'dp 


' du x du y du z 


dx + dy 


+ 


dz 


< 51 / 


dt 


+ pV • u — 0. 


(2,3) 


(2.4) 


Here, we use the chain rule of differentiation to expand the total time derivative into partial 
derivatives with respect to position and time: 

dp _ dp dx dp dy dp dz dp 

dt dt dt dx dt dy + dt dz 


dp 

dt 


+ u • Vp. 


(2.5) 


Since p in (2.5) can represent any scalar function or component of a vector function, the 
convective derivative can be generally defined as 


d 

dt 


d 

dt 


+ u • V. 


(2,6) 


Note that the total time derivative refers to the changes occurring in a fluid element as we 
follow it in its motion, and when it is used explicitly in the equations of motion, such a for- 
mulation is called Lagrangian fluid mechanics. In Lagrangian fluid mechanics, an individual 
velocity is attached to each fluid element as it moves about. Alternatively, if (2.6) is used 
to convert all total time-derivatives to partial differential form, then we have Eulerian fluid 
mechanics, in which we focus on the behavior of a fluid as it moves past fixed points of space. 
In Eulerian fluid mechanics, velocity is a field that has values identified with points in space, 
rather than to individual particles of fluid. In Lagrangian fluid mechanics, one may consider 
the small volume SV as being attached to the fluid element, while in the Eulerian view, 
the small volume 8V is associated with a fixed location in space, and fluid elements move 


2.2. THE NAVIER-STOKES EQUATION 


7 


through it. It is this latter viewpoint that is most common and will be generally adopted 
here (although the Lagrangian viewpoint has its occasional uses). 

Finally, after combining (2.4) (omitting the non-zero factor SV) with (2.5), we get the 
continuity equation: 

^ + V-(pu) = 0. (2.7) 

This is the first of the basic equations of fluid mechanics, and defines how density is related 
to velocity in a compressible fluid. 

Here, however, we wish to consider incompressible flow. In this case, the density p = 
constant, and (2.7) becomes 


Vu = 0. (2.8) 

If (2.8) holds, the velocity field u is often said to be solenoidal, drawing on magnetic termi- 
nology. 


2.2 The Navier-Stokes Equation 

The next basic equation determines the evolution of fluid momentum density. The momen- 
tum of an infinitesimal fluid element, of mass 5m — pSV = constant, is p = 5mu. Applying 
Newton’s second law, the time rate of change of p is due to whatever forces F are imposed: 
dp/dt = F. Upon using (2.6), this becomes 

P 5V + u • Vu) = F. (2.9) 

Now, the various forces that make up F, and thereby affect fluid motion, must be determined. 

First, there is the pressure p in the fluid, which is a force per unit area applied to 
the six faces of the cube of volume SV = 5x5y5z , and the difference in pressure between 
opposing faces will give the net force due to pressure in the direction normal to those faces. 
For example, let the two faces with area 6y5z and normal along the ^-direction have x- 
coordinates x and x + 5x. The component of force in the ^-direction due to pressure points 
from high to low pressure; this has the limiting form: 

F p , x = flm -\p(x + 6x) -p(x)}5y5z 

= 5V. (2.10) 

ox 

The components of the pressure force along the y- and z-directions, F PtV and F p<z , respectively, 
can be determined in a similar manner. Thus, the total force on 5V due to pressure is 

F = 

r p 


— Vp SV. 


( 2 . 11 ) 
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The pressure p is an example of a stress on the surfaces of a fluid element, whose (negative) 
gradient produces a net force. 

Another force may arise due to shear stresses, which are also forces per unit area and 
which occur when neighboring fluid elements move at different velocities. This may be 
thought of as a frictional force between the surfaces of adjoining faces of fluid volumes. 
Consider a cubical fluid volume 5V centered at x = ( x , y, z ) and an adjoining one 5Vj of 
the same volume centered at x' = (x + 5x,y,z). The shear stress cr x is proportional to 
the velocity difference between x' and x (the shear) divided by the distance 5x, with a 
proportionality constant p (the dynamic viscosity ): 

/ X u(x') - u(x) d\x 

o’xfri) = lim p v = P a • (2.12) 

<5x->0 dX OX 

The shear stress can be visualized as manifesting itself on the surface between 8 V and 5V\, 
such that a force acts at the position xi = (x + 5x/2,y,z). Similarly, consider another 
volume 8 V 2 located at x" = (x — 5x,y,z)\ by changing +5x to — 5x in (2.12), and defining 
X 2 — (x — 5x/ 2, y, z), we have 


/ \ ,. u(x") — u(x) chi 

<r x (x 2 ) = bm p = — p . (2.13) 

V 6 x^ 0 ^ Sx ^ dx v ’ 

The net shear force F S)a; acting on the x-directed faces (of area 5y8z ) of the fluid volume 
8V, is 

F s ,*(x) = 


Although p may, in general, be a function of position, here we will assume p = constant. In 
this case, 


lim [<t*(x x ) + <r x (x 2 )] Sy8z 

ox-8-0 


d<j 2 
dx 

L (Q~J 6V ■ 


SxSy8z 
' du 


(2.14) 


F„,!x) = SV. (2.15) 

Since the procedure for finding F SiX can similarly be applied to find the shear forces F SiV and 
F s>z due to viscous stresses acting on the faces of 5V normal to the y- and ^-directions, we 
see immediately from (2.15) that x — > y produces F s>y and x -4 z produces F Sj2 . Adding all 
these shear forces together yields the total force F s due to viscous shear on a fluid element 
5m: 


F s — F S)X + F SiV + F Sj2 — pV 2 u <5V. 


(2.16) 
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Thus we have found the second of two forces that act on the surfaces of an incompressible 
fluid element, and thereby affect its motion. 

Putting (2.11) and (2.16) into (2.9), we arrive at the Navier-Stokes equation : 

= — Vp + pV 2 u + f . (2.17) 

Here f denotes any additional force densities related to body forces F = fSV that may 
be acting on a fluid element. In classical physics, the only available candidates for F are 
gravitational or electromagnetic forces. Gravity is only important when compressibility is, 
so it is not of concern here. The electromagnetic force is important in a magnetofluid, and is 
due to the presence of electrical currents and magnetic fields. It will be discussed presently. 

The pressure p appears to be a new independent fluid variable in (2.17), and would require 
another time-evolution equation to describe it, if the fluid were compressible. However, for 
incompressible flow, p can be determined by taking the divergence of (2.17), using (2.8) to 
produce 


'du 
. dt 


+ u • Vu 


V 2 p = V • (f — pu • Vu) . (2.18) 

Thus, p is found for incompressible flow by solving the Poisson equation (2.18), which depends 
only on the instantaneous values of f and u (since p is constant). 

As Helmholz’s theorem tells us, a vector field is determined by its divergence, curl and 
boundary values, and we can apply this to further define the equations for u. Here, periodic 
boundary conditions will be in effect, and the divergence of u is zero, by (2.8), so what 
remains is to find an equation for the curl of u, which is called the vorticity: w = V x u. 
Taking the curl of (2.17), and using some well known vector identities, produces the vorticity 
equation: 

= V x (u x w + + vY7 2 oj. (2.19) 

dt V PJ 

Here, the kinematic viscosity is u — p/p and is also a constant. Electromagnetic forces are 
represented through the introduction of the appropriate f in (2.19), which is our next topic. 


2.3 Magnetohydrodynamics 

When a fluid is able to conduct electricity, the subject of fluid mechanics expands into what is 
known as MHD or magnetofluid dynamics. Since electromagnetism is now included, perhaps 
a good place to start a discussion is with Maxwell’s equations, the governing equations of all 
electric fields E and magnetic fields B. In SI units, and in free space, these equations are 


V B 


0 


( 2 . 20 ) 
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V -E 


V x B 

V x E 


Pe 


e 0 


Poj + 


i dE 

c 2 dt 


dB 

dt 


( 2 . 21 ) 

(2.22) 

(2.23) 


The constants e 0 and p, 0 are the electric permittivity and magnetic permeability, respectively, 
of free space, while c is the speed of light, p e is the electric charge density and j is the electric 
current density. The last term on the right side of (2.22), divided by p 0 , is called the 
displacement current. 

For incompressible fluid mechanics, characteristic flow velocities U 0 are small with respect 
to the speed of sound, and thus very small with respect to the speed of light c. If L 0 
is a characteristic length, and T 0 = L 0 /U 0 is a characteristic time, then (2.23) tells us 
that nominal magnitudes of E and B (call them E 0 and B 0 . respectively) are related by 
B 0 U 0 . Using this, the relative value of the last term on the right side of (2.22) is 


\c~ 2 dE/dt\ U 2 
|V x B| ~ c 2 


(2.24) 


Thus, the term c 2 <9E /dt is negligible compared to the other terms in (2.22). The equation 
(2.22) then reduces to a defining relation between current and field: 


p 0 j = VxB. (2.25) 

[As a historical aside, Maxwell generalized (2.22) from the known time-independent rela- 
tionship (2.25), so that taking the divergence of (2.22) and using (2.21) would lead to a 
continuity equation for electric charge. This discovery of the ‘displacement current’ com- 
pleted the mathematical foundations of electromagnetism, which were thereafter called, in 
his honor, Maxwell’s equations.] 

While (2.24) indicates that the time rate of change of E is negligible in MHD, this is not 
necessarily true for E itself. In fact, we can find an expression for E, and through this, an 
expression for the charge density p e . First, we can introduce a vector potential A because of 
( 2 . 20 ): 


V • B = 0 -> B = V x A. 


(2.26) 


Here we may choose a condition (called a gauge condition) on the divergence of A. We 
choose 

V-A = 0. (2.27) 


Second, we use (2.23) and (2.26) to introduce a scalar potential (p: 


V x E = 


d B 

dt 


-> Vx E + 


d A' 

dt 


= 0 


^ E+ at = ~ Vv ' 


(2.28) 
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The electric field E is therefore determined by <p and A: 

f) 

E = - m . (2.29) 

Finally, taking the divergence of (2.29), then using (2.21) and (2.27), yields 

V-E = -VV = Pe • (2-30) 

To reiterate, although <9E /dt is neglected in the MHD approximation, E and p e are generally 
nonzero. 

Now we determine the form of E. Let the fields which are seen in the Lagrangian reference 
frame attached to a moving fluid element be denoted by E' and B', and the (Eulerian) 
fields as seen at fixed points of space by E and B. These fields are connected by Lorentz 
transformations which, for the low velocities of interest here, become 

E' = E + uxB and B' = B. (2.31) 


In addition, the current j in the reference frame attached to a moving fluid element can be 
given by the simple Ohm’s law: j = crE', where a is the electrical conductivity of the fluid 
(which will be assumed constant here). Upon using (2.31), Ohm’s law becomes 

j = cr(E + uxB). (2.32) 

Solving (2.32) for E and using (2.25) for j produces 

E = r]V x B — u x B. (2.33) 


The constant rj is the resistivity : rj — (yw 0 <x) 1 ■ Combining (2.33) and (2.29) gives the 
evolution equation for the vector potential: 


dA 

dt 


— V</? + u x B — 77 V x B. 


(2.34) 


Taking the divergence of (2.34), and using (2.27) and (2.30), gives the MHD expression for 
charge density: 

pe = -VV = -V-(uxB). (2.35) 


Note that (2.35) and (2.18) are similar in the way they define electrostatic potential tp and 
static pressure p, respectively. 

Putting (2.33) into (2.23), or equivalently, taking the curl of (2.34), gives the magnetic 
field evolution equation : 


dB 

dt 


= V x (u x B) + rjV 2 B. 


(2.36) 
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Note that this equation for the magnetic field B has a similar structure to (2.19). Since 
B is the evolving quantity, and j is defined by (2.25), it is important to recognize that, in 
the MHD approximation, B determines j, rather than the reverse, which occurs in the more 
general case when (2.22) must be used instead of (2.25). 

The electromagnetic force density f appearing in (2.19) is 

f = j x B. (2.37) 

Here, this will be called the magnetic force density. It is essentially the vector sum of the 
Lorentz forces qv x B acting on the individual particles of charge q and velocity v which 
collectively constitute the ‘fluid’ in SV. Placing (2.37) into (2.19) gives the vorticity equation : 

d “ = Vx(uxw+ 1 jxBl+ (2.38) 

dt V P J 

There are two dissipation coefficients in (2.38) and (2.36): the viscosity v and the resistivity 
7], If no energy is input into the magnetofluid, the presence of u and rj ensures that any 
initial energy will decay towards zero, as will be shown in the next chapter. First, let us look 
at the non-dimensional form of the equations of motion. 


2.4 Non-dimensional Equations 

The formal procedure for producing non-dimensional equations begins with the following 
assignments: 

u = U 0 u', B = B 0 B', x = L 0 x', t=fy- (2.39) 

Furthermore let B 0 = y/p 0 pU 0 , which equates characteristic velocity U Q with the so-called 
Alfven velocity, Ua = B 0 /y/p 0 p. Placing (2.39) along with p = constant into (2.38) and 
(2.36) yields 

^ = V'xfu'xw'+j'xB'j + ^V'V (2.40) 

= V' x (u' xB')+ V' 2 B'. (2.41) 

Ot ±i]\/[ 

The Reynolds number R e , magnetic Reynolds number R M , and V' are 

R e = UoL °, R m = U ° L ° V'= a ,. (2.42) 

v r] ax' 

The dimensionless numbers R e and Rm characterize the magneto-flow, and any two systems, 
independent of relative physical size, with the same values of R e and Rm are said to have 
similar flow. 
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We will always work with the equations in non-dimensional form. In this case, it is 
convenient to remove all the primes in (2.40) and (2.41), and also to set R~ l — v and 
Rm = rj , for brevity, where u and rj are now dimensionless numbers. The result is a set of 
non-dimensional equations for the time development of the fields u and B: 


duj 

dt 


Vx(uxw+jxB) + vV 2 u> 


(2.43) 


dB 

dt 


V x (u x B) + t?V 2 B. 


(2.44) 


In addition to these, we also have the non-dimensional relations: 


V • u = 0, oj — V x u 


V • B = 0, j = V x B. (2.45) 

The equations (2.43), (2.44), and (2.45) define the basic equations needed here. If we set 
rj = v = 0, we obtain the equations of ideal magnetofluid turbulence. (The equations of 
simple fluid turbulence arise if we set B — 0.) 

Finally, boundary conditions remain to be specified. In general, if solid surfaces S with 
local normal n are present, then n • u|s = 0 for inviscid flow, while u|s = 0 for viscous 
flow. The magnetic boundary conditions depend on the conductivity and permeability of any 
surface, and also require specification. However, we need not go into boundary conditions any 
further, since we assume here and henceforth that any flow under consideration occurs in a 
periodic box. Periodic boundary conditions will be satisfied using spatial Fourier expansions 
of u and B, as will be seen in the next chapter. 


2.5 References for Further Reading 

Perhaps the best all-round book on fluids (with separate chapters on ideal fluids and turbu- 
lence) is [Landau 87]. 

MHD is discussed by the same authors in [Landau 84, Chap. VIII]. 

A derivation of the MHD equations from the underlying plasma kinetic equation is found in 
[Nicholson 83]. 

Boundary conditions on electromagnetic fields are discussed in [Jackson 75]. 

Suggestions for references on the subject of turbulence have been given at the end of the 
previous chapter. 
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There are many more books that deal with various aspects of fluids, magnetofluids, and 
turbulence, as a library or net search will show. Some of these will be mentioned in the 
chapters ahead. 

Finally, a good mathematical reference, on Helmholz’s theorem, among other topics, is 
[Arfken 00]. 


Chapter 3 

Fourier Representation 


Here, it is assumed that the basic character of a given turbulent flow is invariant with respect 
to spatial translation; such flow is termed homogeneous. Our principle concern is with the 
behavior of the flow within an arbitrary region, far removed from any boundaries, somewhere 
in a large extent of fluid. In particular, we assume the large region can be divided into boxes 
and is periodic from one box to any adjoining ones. Alternatively, one can consider the 
flow as occurring in an unbounded but finite space called a 2-torus, for 2-dimensional (2- 
D) flow, or a 3-torus, for 3-dimensional (3-D) flow. In either case we have a finite area or 
volume of space whose points are defined by position vectors x, where x = a;x + yy for two 
dimensions, and x = xx + yy + zz for three dimensions (x, y, z are unit vectors in the x, 
y, z directions, respectively). The components of x take on values mod(27r), so that their 
domains are 0<rr<27r,0<y<27r, and 0 < z < 27T, which are suitable for examining 
spatially periodic solutions of the non-dimensional equations of motion. 

For brevity, we will set X\ — x, = y, and £3 = z, as well as ei = x, §2 = y, and 63 = z. 
Then we can write 
D 

X = Y^XiCi- (3-1) 

i — 1 

Here, D = 2 or 3, signifying a vector in a 2-D or 3-D space, respectively. Also, 0 = (0, 0) or 
(0,0,0), respectively, in 2-D or 3-D. 

3.1 Discrete Fourier Transforms 

Since there is periodicity in each spatial direction, we can represent the magnetofluid vari- 
ables u> and B, assuming sufficient smoothness, in terms of discrete Fourier series. We begin 
by defining a finite set K$ of wave- vectors k (where / is the set of all integers): 

D 

k = kj ki e I, D = 2 or 3 

i= 1 

I<d = {k | -N/2 < kj < N/2, j = (3.2) 
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(Typically, N will be a power of 2: N = 2 m .) 

Next, we define a finite set Xp of position vectors x with discrete coefficients Xi = 
27 T rii/N: 


27T 

x = N 

D 

£ ni ej, ni el, 
2 — 1 

D = 2 or 3 


Xg = {x 

1 0 < rij < N — 1, 

j = l, 

(3.3) 


The sets Xp and Kp enable the definition of finite, discrete Fourier transformations: 

«(x) = £ w(k) e ik x , x € Xg (3.4) 

keK% 

B(x) = £ B(k) e ikx , x € Xp. (3.5) 

k eK% 

The inverse transformations are 


(i(k) = 

£ «(*) e~ ikx , 

xexg 

k e i< n d 

(3.6) 

B(k) = 

L E B(x) e"’ k x , 

X6-XJ? 

k e Kg. 

(3.7) 


Relations (3.4) and (3.5) are said to transform the variables from k-space to x-space, and 
(3.6) and (3.7) reverse this transformation. 

In the above, it will be noticed that the number of discrete x-space points x e Xp is 
N D , while the number of discrete k-space points k e Kp appears to be (N + 1) D . However, 
from (3.6) and (3.7) (and suppressing explicit time-dependence for brevity), it is clear that 

u>(k* + ) = c5(kL) and B(k* + ) = B(k!_), i = (3.8) 

where the k± are defined for D = 3 by 


N 

k± = ± 2 ®1 + & 2©2 + ^363 

(3.9) 

,2 , . 

k ± = fciei ± 2 e 2 + k 3 e 3 

(3.10) 

3 N 

k ± = kiei + k 2 e 2 ± e 3 . 

(3.11) 

For D — 2, we use only (3.9) and (3.10), with k 3 = 0. 


Thus, due to (3.8) we have equivalent points k!,. ~ k*_, i = 1, . . 

. , D, in the set Kp , so 


that the number of independent points in k-space is decreased by one in each dimension. 
The total number of points needed for invertible Fourier transformation is therefore N D in 
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both x-space and k-space. Note that (3.8) explicitly shows that the k-space representations 
of the magnetofluid variables, as well as those in x-space, are periodic. 

Also, it is important to note that values of x other than those defined by (3.3) can be used 
in the Fourier transformations. Shifting all the x € X% by the same amount x x+d merely 
multiplies the complex Fourier coefficient associated with k by a factor a(k) = exp(ik • d), 
thereby changing the phase of the coefficient. Conversely, multiplying coefficients by a(k) 
allows the values of x not in Xp to be determined. Thus, a truncated discrete Fourier 
transformation is actually a continuum representation of a physical variable, although only 
to the level of resolution allowed by the discrete number of k 6 Kp. The x € X$ may 
be viewed as a set of sampling points of the continuum, sufficient to reproduce the values 
at any point in the periodic space. The continuum nature of the representation also means 
that spatial derivatives can be evaluated exactly: V /(x) —> ikf(k). This makes numerical 
methods based on Fourier representations (whenever they can be applied) inherently more 
accurate than finite difference methods, which have no underlying continuum nature. 

In k-space, (2.45) becomes 

k • u(k) — 0, a? (k) = ikx u(k) 

k-B(k) = 0, j(k) = ikxB(k). (3.12) 

Also, from the above, we may write 

u(k) = ik~ 2 k x w(k) and B(k) = ik~ 2 k x j(k). (3.13) 

Obviously, k • u?(k) =0 and k • j(k) = 0, because of (3.12). 

Next, we consider information content in the real- valued variables u(x) and B(x), and 
complex- valued variables u(k) and B(k). In 3-D x-space, there appear to be 6N 3 real values 
(the 3 components of u and the 3 components of B at each discrete value of x), while in 
3-D k-space, there appear to be 12 N 3 real values (the 3 complex components of u and the 3 
complex components of B at each discrete value of k, taking (3.8) into account). Similarly, in 
2-D x-space, there appear to be 4 N 2 real values (the 2 components of u and the 2 components 
of B at each discrete value of x), while in 2-D k-space, there appear to be 8 N 2 real values 
(the 2 complex components of u and the 2 complex components of B at each discrete value 
of k, again taking (3.8) into account). However, since u(x) [or equivalently, w(x)] and B(x) 
are real-valued, it is evident from (3.6) and (3.7) that 

w*(k)=w(-k) and B*(k) = B(-k). (3.14) 

This reality condition reduces the number of components (real and imaginary) which may 
be independent by half in 3-D k-space, to 6 N 3 , which is commensurate with the 6 N 3 real 
components in 3-D x-space. Similarly, the number of components (real and imaginary) in 
2-D k-space is also reduced by half, so that only 4iV 2 are independent, which is commen- 
surate with the AN 2 real components in 2-D x-space. In fact, the relations (3.14) are used 
to minimize computer storage in subroutines which perform Fourier transformations of real 
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variables, so that (3.14) are (almost) explicitly enforced in any computer simulation of ho- 
mogeneous turbulence that uses Fourier transforms. (Almost, because for a small number 
of k, coefficients for — k are also utilized in Fast Fourier Transform, or FFT, subroutines.) 

There are further reductions possible in information content, which are easier to see in 
k-space, than in x-space. First, a reduction is brought about bv the solenoidal relations 
k • u(k) = 0 and k • B(k) = 0 in (3.12). These solenoidality conditions reduce the number 
of dynamically independent components of u.(k) and B(k) by an additional factor of 2/3 for 
D = 3 and by a factor of 1/2 for D = 2. The 2-D vectors u(k) and B(k) can be represented 
by exactly one complex function each, thereby explicitly incorporating the 1/2 reduction, 
while in 3-D, the vectors u(k) and B(k) generally have three complex components each. 
However, these three complex components are kimematically linked by (3.12), so that of the 
three, only two are independent. 

Using these observations, we see that in 3-D the number of independent real and imagi- 
nary parts of either u(k) or B(k) is 2N Z rather than 3 iV 3 , for the N z independent k G , 
while for 2-D, this number is N 2 , equal to the N 2 independent k G KI? . However, we do not 
need to use all of the possible k G Kp ; we can instead algorithmically include only those ±k 
from a chosen subset S C Kp . The structure of possible subsets S C K will be further 
discussed after an examination of the form of the dynamical equations in k-space. 


3.2 Evolution Equations in k-Space 

3.2.1 Three-Dimensional Equations 

The evolution equations of vorticity and magnetic field in k-space are found by placing (3.4) 
and (3.5) into (2.43) and (2.44). The result is, using the independence of the e lk x , 

= S(u, w; k) + S(j, B; k) — vk 2 u:(k) (3.15) 

= S(u, B; k) — 77& 2 B(k). (3.16) 

Here, the nonlinear terms denoted by S are vector convolutions: 

p+q=k 

S(u, B; k) = ikx £ 

k, p.q eS 

The double summation Ek^qes (3.17) is over all wave-vectors p and q that satisfy 
p + q = k, where ±k, ±p, ±q G S C Kp. 

Thus, Fourier coefficients, such as u(k), are associated only with a finite number of 
k G S C K$, and not with all k G K Notice that, using (3.17), the right sides of (3.15) 
and (3.16) are identically zero for k = 0. Thus, we may set u.(0) = 0 by Galilean invariance 
(a value it holds for all t). while the constant value of B(0) may be chosen to be zero or 


[u(p) x B(q)] . (3.17) 


dui( k) 

dt 

c/B(k) 

dt 
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nonzero to reflect the absence or presence of an externally imposed mean magnetic field. In 
order to ensure homogeneity, we set a>(0) =0 to reflect no overall rotation in the fluid. Also, 
since k- S(u, B; k) = 0, then, if (3.12) holds initially, (3.15) and (3.16) show that it will hold 
for t > 0. 

The right side of (3.17) contains a vector triple-product that can be expanded as 

p+ q =k 

S(u,B;k) = i ^ k x [u(p) x B(q) 

k, p,q eS 

p+q=k 

= * £ [p • B(q) u(p) - ( 

k, p,q eS 


•u(p)B(q)]. (3.18) 


The last step follows from k = p + q and the solenoidality (3.12) of u and B. If p and 
q are collinear, then, for example, q • u(p) ~ p • u(p) = 0, and these wave vectors do not 
contribute to the right side of (3.18). Thus, in order to have a nonlinear interaction between 
modes with wave vectors p, q G <5, we must have pxq/0. 

The 3-D vector convolution sums S defined in (3.17) can also be written as 

S(u, B; k) - C(k) • Q(u, B;k), (3.19) 


where the dyadic C(k) and the vector Q are defined by 

e* • C(k) • e.j = Qj(k) = ie imj k m . (3.20) 

p +q= k 

Q(u,B;k) = £ u(p) x B(q). (3.21) 

k >p,qe-s 

Here, the summation convention is used: Ci m jk m = timjkmi so that repeated indices 

imply a summation. The tij m are the components of the completely alternating tensor , or 
Levi-Civita symbol, which is defined by 


1, if (i,j,k) are cyclic permutations of (1,2,3) 

= —1, if (i,j,k) are anticyclic permutations of (1,2,3) 

0, otherwise. 


( 3 . 22 ) 


In addition to the curl operator C(k), there are several other important dyadics. The 
first of these will be denoted by P(k). It arises when we use (3.13) together with (3.12), to 
get 

u(k) = ik~ 2 k x u>(k) 

= — &T 2 k x [k x u(k)] 

= —AT 2 kk ■ u(k) — /c 2 u(k)| 

= P(k) • u(k). 


(3.23) 
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This projection operator P(k) is the dyadic defined by 

P(k) -'“F = 0-24) 

Here, the second important dyadic is I, which is the identity operator. Next, we define the 
adjoint curl (or ‘uncurl’) operator C(k): 

C(k) = k~ 2 C(k). (3.25) 

The projection, curl and adjoint curl operators satisfy 

k-P(k) = P(k)-k 
k • C(k) = C(k) • k 
k • C(k) = C(k) • k 
P(k) • P(k) = P(k) 

C(k)-P(k) = P(k) • C(k) 

C(k) - P(k) = P(k) • C(k) 

C(k) • C(k) = C(k) • C(k) 

C(k) • C(k) = fc 2 P(k) 

C(k) • C(k) = AT 2 P(k). 

Note that C(k) is not quite the inverse of C(k). 

We can use the operators just defined, to rewrite (3.13) as 


u(k) = C(k) • u?(k) = — u;(k) • C(k) (3.27) 

B(k) = C(k) • j(k) = — j(k) • C(k) (3.28) 

and also to rewrite (3.15) and (3.16) as 

d ^dt^ = + ^ B k ^l ~~ ^ 2ci, ( k ) ( 3 - 29 ) 

dB ^ = C(k)-Q(u,B;k)-77A: 2 B(k). (3.30) 

Now, we can apply C(k) to (3.28) to arrive at a definition of the k-space stream vector tj> 
and magnetic vector potential A: 

t/>(k) = C(k) • u(k) = k~ 2 u>(k) (3.31) 

A(k) = C(k) • B(k) = AT 2 j(k). (3.32) 



3.2. EVOLUTION EQUATIONS IN K-SPACE 


21 


Using (3.32 and (3.28), along with (3.21) and (3.26), and applying C(k) to (3.29) and (3.30) 
gives 


dU ^ = P(k).[Q(u,w;k) + Q(j,B;k)]-i/fc 2 u(k) (3.33) 

dk }^ = P(k) ■ Q(u, B; k) — r/k 2 A(k). (3.34) 

at 

Equations for ijj and j could be similarly generated, if desired. 

It has already been mentioned that the magnetic field may contain a nonzero component 
B(k) for k = 0. Let the constant mean magnetic field be B 0 = B(0), and let b(k) = B(k) 
for k ^ 0, with b(0) = 0, so that (3.21) becomes 

Q(u, B;k) = Q(u, b;k) + u(k) x B 0 . (3.35) 

Also, we will define a(k) = ik~ 2 k x b(k) = &r 2 j(k). Since j(0) = 0, then Q(j,B;k) can 
be obtained by letting u -» j in (3.35). These results allow' the evolution equations (3.29), 
(3.30), (3.33), and (3.34) to be written 


dUJ }t^ = C ( k ) ' + Q(j>b;k)] 

+ i(k • B 0 )i(k) - uk 2 u)( k) 


db(k) 

dt 

du(k) 

dt 


da(k) 

dt 


C(k) • Q(u, b; k) + i(k • B 0 ) u(k) - r?fc 2 b(k). 

P(k) • [Q(u,w;k) + Q(j,b;k)] 

+ i ( k • B 0 ) b(k) — uk 2 u(k) 

P(k) • Q(u, b; k) + i(k • B 0 ) T/>(k) - rjk 2 a{k). 


(3.36) 


(3.37) 

(3.38) 

(3.39) 


Many more equations can be generated through the application of C(k) and C(k) to these 
equations. However, for our purposes here, the above will be sufficient. 


3.2.2 Two-Dimensional Equations 

The equations (3.36), (3.37), (3.38), and (3.39) are fully three-dimensional but can be reduced 
to a 2-D form as follows. In 2-D x-space, (2.45) allows us to write vorticity and current as 

w(x) = w(x)z and j(x)=j(x)z, 


(3.40) 
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where x = xx + yy. In 2-D k-space, using (3.12), (3.31), and (3.32) 
relation between Fourier coefficients: 

leads to the following 

Y 

ii 

= k?ip(k)i 

— k^ipifk) 

(3.41) 

C— .i 

II 

^.1 

= k 2 a( k)z 

— /j 2 a(k) 

(3.42) 


where k = k x x + k y y, ip is the scalar stream function and a is the scalar magnetic potential. 
In the 2-D case, (3.13), along with (3.41) and (3.42), leads to the following expressions for 
the velocity and magnetic fields in k-space: 


u(k) = i k x z ip(k) and b(k) =ikxza(k). (3.43) 

Thus there is exactly only one independent complex component in u(k) and one in b(k), 
a reflection of the fact that they are 2-D solenoidal vector fields. A 2-D magnetofluid can 
thus be described by ip( k) and a(k), whose evolution equations follow from equations (3.43), 
(3.36), and (3.39): 

p+q=k r - _ , 

i Y Z'Pxq[^(p)w(q)+i(p)a(q)] 

k,p,q£5 

+ i(k • B 0 ) j(k) — uk 2 Cj(k) (3.44) 


d<u)(k) 

dt 


p+q=k 

-i Y z • p x q V’(p) fi(q) 

k,p,q6<S 

+ i(k • B 0 ) ip(k) — rjk^aik). (3.45) 

Here, we have used u)(k) = k 2 ip( k) and j(k) = k 2 a( k), as given in (3.41) and (3.42). Since p 
and q are 2-D vectors, pxq = ( p x q y — p y q x )z and z • p x q = p x q y — p y q x . Also, remember 
that the double sum in (3.45) is over all p and q such that p + q = k and k, p, q € S C K? . 
As in the 3-D case, equations (3.44) and (3.45) clearly indicate that only wave vectors p and 
q such that p x q/0 can interact nonlinearly and thereby affect the time-evolution of the 
Fourier mode k. 


da( k) 

dt 


3.2.3 Alfven Waves 

Here we consider the dynamical equations when B 0 is relatively large. Let us define the 
Alfven frequency for a mode k as 

C(k) = k B 0 . (3.46) 

In the 3-D case, if all the terms on the right sides of equations (3.37) and (3.38) except the 
terms containing £(k) are dropped, these equations become 

db jf = *C(k)u(k) 


(3.47) 
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(3.48) 


(3.49) 


(3.50) 

Using (3.47) produces the corresponding velocity coefficient solutions: 

u (±) (k) = ±b (±) (k). (3.51) 

The x-space waves corresponding to these Fourier coefficients are 

b^(x,t) = ±u^\x,t) = b£^(k) exp[ik • (x ± B 0 £)]. 

(3.52) 

The linear waves b^ (x, t ) = (x, t ) are called Alfven waves. Furthermore, the waves 

b (~)(x,t) and u ( _ )(x,£) propagate in the direction of the corresponding k, with speed B 0 , 
while b( + )(x, t) and u^ + ^(x, t) propagate in the opposite direction at the same speed. In the 
linear approximation, Alfven waves do not interact, but if B 0 is no longer ‘large,’ then they 
do interact, in which case Alfven waves may lose their individual identities. 

The linear analysis above produces pairs of solutions b (±) (x, t) and u (±) (x, t), which 
are coplanar. Thus, these results go over directly into the 2-D case, where the equations 
corresponding to (3.47) and (3.48) come from the linearization of (3.44) and (3.45), with 
v = r/ = 0. Using (3.43), (3.50), and (3.51), along with and (3.41), yields the ideal 2-D MHD 
Alfven wave solutions in k-space: 

a^(k) = ±-0^(k) = a^(k) exp[±*£(k)t]. (3.53) 

The corresponding x-space equations are 

a^(x,t) = ±?/)^(x, t) = a^(k) exp[*k • (x ± B 0 t)]. 

(3.54) 

Although, in general, Alfven waves may lose their individual identities in the presence of 
nonlinearity, some remnants of Alfven wave-like behavior may persist in MHD turbulence. 
In fact, this will be apparent when numerical results are presented in a later chapter. 


du(k) 

dt 


*C(k)b(k). 


These can be combined into one linear, second-order differential equation: 

d 2 b(k) 


dt 2 


= -< 2 (k)b(k), 


which has two solutions: 


b (±) (k) = b^, ±) (k) exp[±*C(k)i]. 
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3.3 Interacting Fourier Modes 

In the previous section, nonlinear evolution equations for the Fourier modes u>(k) and b(k) 
were developed, for both the 2-D and 3-D cases. Once we choose a grid size for a numerical 
simulation, we have a finite number of k e Kp to choose from; in fact, we choose ±k € 
S C Kp, and the choice of S is essentially arbitrary. For example we could choose to work 
with only one k = kx ^ 0 (and — k, of course), setting the Fourier modes associated with 
all other k £ S C Kp equal to zero. In this case, only the mode associated with ±ki is 
time advanced, but there is no nonlinear evolution, since all other modes are frozen with 
zero values. This is, of course, a trivial case, but there are many other possible subsets S of 
Kp that can be chosen. The simplest nonlinearly interacting system is one associated with 
a subset of Kp containing three wave vectors (and their negatives). All other interacting 
subsets of Kp will be unions of these basic interacting triads. 

However, choosing one or an arbitrary union of interacting triads may impose an un- 
wanted anisotropy on the model system, if what is really desired is to have a set of wave 
vectors which give no intrinsic directionality to the model system. Restricting the choice of 
a subset S of Kp so that no particular direction in space is favored is called an isotropic 
truncation. In what follows, we specifically discuss these two choices: interacting triads and 
isotropic truncation. 

3.3.1 Interacting Triads 

The basic nonlinear interaction in ideal turbulence is quadratic, as was seen explicitly in 
the last section. In terms of nonzero wave vectors k, p,q £ Kp , D = 2,3, Fourier modes 
associated with nonzero wave vectors ±k, ±p, ±q can nonlinearly interact only if they form 
a fundamental interacting triad: k = p + q, where p x q / 0. These are the minimal 
interacting subsets of the Kp (any other interacting subset is a union of these fundamental 
subsets). 

An interacting triad T C Kp thus consists of three independent wave vectors k, p, q € 

T< n - 
1s -d . 

T = {±k,±p,±q} 

k = p + q, pxq/0. (3.55) 

Note that, because of (3.14), we must have both a wave vector and its negative in T (however, 
remember that k and — k identify the same mode). Also, note that although T has wave 
vectors that satisfy k = p + q, it does not contain k' = p — q. If a model system were to 
include only modes associated with wave vectors in T, then only these would advance in 
time and all other modes would remain zero, having been set so initially. (In atmospheric 
science, truncation of a buoyant fluid model to three modes led to the well-known Lorenz 
attractor.) 

Using the definition, (3.55), consider the two sets of interacting triads T\ and T 2 '. 

T\ - {±ki, ±pi, ±qi} 
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T 2 = {±k 2 ,±P2,±q2}- (3-56) 

If 71 D T 2 = 0, then the subsets Ti and % do not interact with one another. However, 
assuming that at least one independent wave vector in 7} is equal to one in T 2 (for example, 
pj = p 2 ), then the subset S = T\ CT 2 is also a set of interacting wave vectors which contains 
T\ and T 2 : 

S = {±ki,±Pi,±qi | i = 1,2} (3.57) 

For example, consider the following interacting triads of 2-D wave vectors: 

T\ = {±(1,0), ±(0,1), ±(1,1)} 

r 2 = {±( 1 , 0 ), ±(o,i), ±(i,-i)> 

r 3 - {±( 2 , 1 ), ±(i,i), ±(i,o)} 

r 4 = {±( 2 , 0 ), ±(0,2), ±(2,2)} (3.58) 

The first three triads, 7}, i = 1,2,3, share common elements and pair-wise unions of any two 
of them will yield larger subsets of interacting wave vectors. However, T^OTi — 0, i = 1,2,3, 
so that any union Si = Ta U 7} , i — 1,2,3 will consist of two non-interacting triads, rather 
than a fully interacting subset. In this last case, Ta and any of the 7 1, i = 1,2,3, though 
formally united, evolve independently of one another. 

We can thus use the subsets 7} C K% (with N > 2) of interacting triads to define larger 
sets, (Si and S 2 , such that <Si D S 2 = 0. These two sets Si and S 2 do not interact with one 
another, and if there are no other interacting sets allowed in the numerical simulation which 
would link them, then they will evolve independent of one another. For N > 2, any nonzero 
k G Kp can always be written as k = p ± q, with nonzero p, q G K%, where p and q satisfy 
P x q ^ O' 

To show this, consider a nonzero, 3-D wave vector: k = k x Sc ± k y y ± k z z. If at least two 
of the ki , i = x.y, z, are definitely nonzero, say k x and k y , then we can immediately write: 

p = k x x, q = ky-y ± k z i -> k = p ± q. (3.59) 

Here we have p 2 , q 2 < k 2 and pxq = —k x k z yi ± k x k y i ^ 0. 

If only one of the k z , i — x,y, z, is nonzero, and the other two are zero, say k x — ±k 
(k > 0) and k y = k z = 0, so that k = ±/cx, then we have two cases: 

p = (k - l)x ± y, q = x - y k = p ± q (k x - k > 0) 

p = ( — k ± l)x - y, q = — x ± y -± k = p ± q (k x = -K< 0). 

(3.60) 
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In both cases ( k x = ±«), we have p 2 = k 2 — 2 (k — 1), q 2 = 2, and pxq = — /cz ^ 0. Here, 
if k > 1 , then p 2 , q 2 < k 2 = k 2 , while if k = 1 , we have p 2 = k 2 — 1 < q 2 = 2. 

3.3.2 Isotropic Truncation 

If we were to pick an upper limit to the magnitude of a wave vector, say k max , where 
kmax > 2, then the results (3.59) and (3.60) indicate that any vector k that satisfies k 2 < 
kmax < (2V/2) 2 , can be decomposed into two wave vectors p and q, such that k = p + q, 
with p 2 ,q 2 < k ^ nax . This gives us another type of subset, /C C K% (N = 2 M > 2), which 
also has all triads linked together. Thus, the isotropic truncation of K$ produces 

K = {k I k e K%, k 2 < k 2 ma , < (N/: 2) 2 }. (3.61) 

In defining K, by (3.61), we have constrained the number of dynamically important k so as 
to avoid imposing any initial anisotropy. We have isotropically truncated the domain of k 
by using only those Fourier components with wave-vectors k lying on or inside the sphere 
defined by k < k max < N/2. All coefficients with k > k max are set to zero and kept there 
during numerical simulation ( i.e ., after every time-step). Only coefficients with k < k max 
are allowed to evolve according to the algorithm used to define the numerical simulation. 

The fact that we set k max < N/2 is a small but important point. Consider the Fourier 
modes corresponding to ki = (N/2, 0,0) (and to k 2 = (0,N/2, 0) and k 3 = (0,0, N/2)). 
Using x as given by (3.3), we have e ik J' x = (— l) n t, j — 1, 2, 3, so that modes determined by 
(3.6) and (3.7) are purely real for k = kj , j — 1, 2, 3. A spatial derivative of such a component 
corresponds to multiplication by ik, which transforms the real part of the component into 
an imaginary one, and the imaginary part (which is zero, in this case) into a real part. Thus, 
the Fourier mode of a derivative associated with one of the k j, which is also purely real, 
is identically zero. Since the evolution of Fourier components is a function of their spatial 
derivatives, any interacting set of wave- vectors should exclude the kj, j — 1,2,3 (as well as 
such wave vectors as k 4 = (N/2, N/2, 0), etc.). 

The number of k within an ‘isotropic’ 2-D circle is approximately equal to M 2 ~ 
and within an ‘isotropic’ 3-D sphere is approximately equal to M 3 47r *4axA More 
precisely, M 2 and M 3 are given by 

M d = £ H(k max -k), D — 2,3, H(a) = J’ 

(3.62) 

Since Mo is the total number of Fourier vector modes actually utilized, and since there is one 
independent component per k in 2-D and two in 3-D, then the total number of independent 
components per physical field (u or B) is, in 2-D, M 2 , and, in 3-D, 2 M 3 . 

Clearly, the total number of independent components in the Fourier modes k £ K% used 
to represent each physical field is reduced by isotropic truncation (3.62). In 2-D k-space, 
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this reduction is to F 2 , and in 3-D k-space, to £ 3 , where 

F 2 = M 2 (3.63) 

Fs = 2M 3 «^„<jAf 3 . (3.64) 

The actual value of Fp depends on the choice of /c max , which is set by computational con- 
straints. However, note that the sums (3.4) and (3.5) can still be taken over all N° modes 
k because those coefficients with k > k ma x have values defined algorithmically as zero. 


3.4 Ideal Integral Invariants 


If v = r] = 0 in the evolution equations, then there exist certain integral invariants which 
are analytical constants of the motion. They are termed ‘integral invariants’ because they 
are integrals over the physical space of various bilinear forms. For example, an integral that 
may be invariant is 

(u,b) = * j u(x) • b(x) d B x, D = 2,3. (3.65) 

If we use Fourier transforms similar to (3.5) for u and b, then (3.65) becomes 


(ii,b) 


(2t t) d / 


£ u(k')e ik ' x 
1 

(2tt) 


k 'etc* 


Y b(k) e ikx 1 


keh-% 


= £ W)-b(k) I f e«^d° X 

k,k>eK” ^ > J 

= Y u(k') • b(k) 5(k + k') 

k,k'€Kp 

= Y fl(-k)-b(k) 

keKg 

= Y u*(k) -b(k), D — 2,3. 

keK% 


d D x 


(3.66) 


Above, the Dirac delta function is 

<5(k + k') = [ e i{ - k>+k) x d D x. (3.67) 

(2tt) v J 

Although we started from the integral (3.65) over a continuous x-space, we could also have 
started with a summation over the discrete x-space (3.3), and defined 

(u, b) = * Y u ( x ) ‘ b(x), £> = 2,3. (3.68) 

xeA-£ 
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Then, using transforms similar to (3.5) again, along with the Kronecker delta 

K* = L E e ,(k ' +k)x . D = 2,3, (3.69) 

X6Xg 

would have led once again to the result in (3.66). Although we could now use the term 
‘summation invariant,’ we will stick to the more standard term of ‘integral invariant.’ 

When v = rj = 0, the integral invariants that are found in the magnetofluid equations 
depend on the specific type of flow that is under consideration. First, the flow may be 3-D 
or 2-D. Then, a flow with B = 0 for all time is called an Euler flow, and with B ^ 0, it is 
called an ideal MHD flow. Finally, an ideal MHD flow may have B 0 = 0 or B 0 ^ 0. Thus, 
there are six different types of flow, and they have different sets of integral invariants. These 
six different cases will now be discussed. 

3.4.1 3-D Invariants 

There are three different cases here: 3-D Euler invariants, 3-D ideal MHD (|B 0 | = B 0 = 0) 
invariants, and 3-D ideal MHD ( B 0 — 1) invariants. [Normalization: In the Euler cases, u 
is usually initialized so that (u, u) = 1, while in the ideal MHD cases, u and b are typically 
initialized so that (u, u) + (b, b) = 2, where (•, •) is defined by (3.65), (3.66) or (3.68).] 


3-D Euler Invariants 


If we categorically set B(x, t) = 0 for all t, then we reduce the number of equations by half. 
In this case, (3.36) and (3.38) become the following: 


G?w(k) 

dt 

du(k) 

dt 


= C(k) • Q(u, ci>; k) — vk 2 u}(k) 
= P(k) -Q(u,w;k) -uk 2 u(k). 


If we take the dot product of (3.71) with u*(k), we get 

d u(k) 


u*(k) 


dt 


— u*(k) • Q(u, w; k) — i4r |u(k)|" 


(3.70) 

(3.71) 


(3.72) 


Now, if we change k -4 — k in (3.72), add the resulting equation to (3.72), and sum over k, 
we get 


E 

k£S 


d|u(k)| s 

dt 


+ 2 v ft 2 |u(k)|^ 


= T(u, u, u>). 


(3.73) 
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Here, using (3.21), the triple-product summation T is defined by 

T(u,u,w) = Y [ fi *( k ) - Q(u,w;k) + u(k) • Q*(u,w;k)] 
k es 

= 2 53 u(k) • Q*(u, w; k) 

kes 

k+p+q=0 

= 2 53 fi ( k ) • fi (p) x "(q) 

k, p,q€<S 
k+p+q=0 

= 53 x fi (p) + fi (p) x a ( k )l • 

k,p,q£5 
= 0 . 


(3.74) 


The triple summation q ls° the summa tion over all wave- vectors k, p, q such that 

k + p + q = 0 and ±k, ±p, ±q e <5 C . We have exercised the option of changing the 
dummy index k — > — k in the second step above, and used the fact that the summation on 
the right side can be made symmetric over k and p (and q also, if necessary). Thus, from 
(3.73) and (3.74) we have, using the equivalence & 2 |u(k)| 2 = |u>(k) | 2 , 


d Ei 

dt 


= —2 uLl. 


Here, kinetic energy Ek and the enstrophy Ll are defined using (3.66): 

Ek = 2 ( u > u ) 

= * (">«)• 


(3.75) 

(3.76) 

(3.77) 


If y — 0, it is clear that Ek is a constant of the motion, i.e., invariant. If v ± 0, then, since 
Ll > 0, Ek always decreases, unless Q = 0, in which case Ek — 0, also. 

Next, if we take the dot product of (3.71) with u>*(k), and take the dot product of (3.70) 
(after changing k — » — k), with u(k), use u(k) • C(k)(— k) = u)(k), add the two equations, 
and sum over k, we get 


£ 

kgs L 


d u(k) • w*(k) 

dt 


+ 2 uk 2 u(k) • w*(k) 


= T(w,u,w). 


(3.78) 


Here, as we did in (3.74), we again partially symmetrize the triple-product sum on the right 
side: 

k+p+q=0 

T(w,u,w) = Y [w(k)xw(q) + w(q)xu(k)]'ii(p) 
k,p,qe<S 

- 0. 


(3.79) 
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Thus, (3.78) becomes 


dH K 

dt 


— 2 V Hy. 


(3.80) 


The kinetic helicity Hr and (giving it a name) the vortical helicity Hy are defined, again 
using (3.66): 


H k = 

2 ( u,tJ ) 

(3.81) 

Hy = 

jKVx w). 

(3.82) 


In the 3-D Euler case, Hk is the second integral invariant. (Neither Hk nor Hy are positive- 
definite.) 

In the present case, there are no other bilinear integral invariants because, looking at 
(3.74) and (3.79), it is clear that in order for the triple-product summation 

k-j-p+q=0 

T(u,v,«) = u(k) • v(p) x u?(q) (3.83) 

k, p, q E«S 

to be equal to zero, at least two out of the three vectors u, v, Co must be equal (up to 
a constant factor). Another way of saying this is to note that, by (3.83), T(u, v,u>) is 
antisymmetric in its arguments: 


T(u, = — T(v, u, w) = — T( u, w,v). (3.84) 

Thus, if u and u> are fixed, there are only two choices for v such that T(u, v, u>) = 0: either 
v = u or v = w. As we have just seen, these choices lead, respectively, to the two invariants 
Ek, defined in (3.76), and Hk, defined in (3.81). 


3-D Ideal MHD Invariants 

In this case, B^O, and we have 3-D ideal MHD flows when v — 77 = 0. If we take the dot 
product of u*(k) with (3.38), then the result of this, added to its complex conjugate, gives 


E 

kes L 


d |u(k)| 2 

dt 


+ 2 i (k • B 0 )u*(k) • b(k) — 2 v |d>(k) | 2 
= T(u,u,w) +T(u,j,b). 


(3.85) 


Similarly, if we take the dot product of b*(k) with (3.37), then the result of this, added to 
its complex conjugate, gives 


E 

ke«s L 


dlb(k)l 2 

dt 


— 2 * (k • B 0 )u*(k) • b(k) — 2 77 |j(k)| 2 
= T(j,u,b). 


(3.86) 
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Adding (3.85) and (3.86), and using (3.84), gives 

dE 

= -2{vD + r]J). (3.87) 

In the above equation, fl is given by (3.77), the total energy E is a sum of kinetic energy 
Ek, given by (3.76), and magnetic energy Em, while J is the mean-square current. We have 

E — Ex + Em (3.88) 

Em = ( (b. h) (3.89) 

J = g (j.J)- (3-90) 

It is clear from (3.87) that E, as given by (3.88), is an integral invariant for u = g = 0. 

Next, take the dot product of b*(k) with (3.38), add this to the dot product of u(k) with 
the complex conjugate of (3.37), and use (3.84) to produce 


E 

kes 


d 

dt 


+ (u + rj) k 2 


u(k)b-(k) 


= 0. 


Defining the cross helicity He as 


H c = 


allows (3.91) to be rewritten as 


dH c 

dt 


1 

2 


(u,b) 


(u + rj)H c . 


(3.91) 


(3.92) 


(3.93) 


The cross helicity is obviously an integral invariant for 3-D ideal MHD, where u = rj = 0. 

Notice that E in (3.88) and He in (3.92) are integral invariants for ideal MHD for both 
B 0 = 0 and B 0 ^ 0. Keeping this in mind, the result of taking the dot product of a*(k) with 
(3.37), and adding this to its complex conjugate, again using (3.84), is 


E 

kes 




a(k)b-(k) 


-ij(k B,)i(k)-r(k). 

keS 


(3.94) 


Defining the magnetic helicity H M as 

H m = 2 ( a > b)> (3.95) 

we see that Hm is an integral invariant in 3-D ideal MHD, but only if B 0 = 0. 

There are no other bilinear invariants for 3-D ideal MHD, and the absence or presence 
of a nonzero mean magnetic field splits 3-D ideal homogeneous MHD turbulence into two 
cases: If B 0 = 0, E , He and Hm are ideal invariants, while if B 0 7^ 0, only E and He are 
ideal invariants. We mention again that E, He, and H M , as well as k • u(k) and k • B(k), 
are zero only to machine accuracy. 
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3.4.2 2-D Invariants 

There are also three different cases here: 2-D Euler invariants, 2-D ideal MHD (|B 0 | = B 0 — 
0) invariants, and 2-D ideal MHD ( B 0 — 1) invariants. 

2-D Euler Invariants 

As in 3-D flow, there is no magnetic field present at all in this case. Let us examine (3.44) 
with j = a = 0: 


dCj( k) 

dt 


p+q=k 

* E z-pxq^(p) U) (q) — vk 2 Cb( k). 

k, p, q 


(3.96) 


If we multiply this equation by <j)* n = k n ip*{ k), add the result to its own complex conjugate, 
and sum over k G S C . we get 


E 

kes L 


d 

dt 


+ 2 v k 2 


k n+2 \^{k)\ 2 




(3.97) 


where 


k+p+q=0 

0((/>,-0,w) = 2 i z P x q^(k)^(p)w(q). (3.98) 

k,p,q£5 

In the summation above, the indices k, p, and q are dummies and can be interchanged freely 
to yield relations similar to (3.84): 


= -6»(^,^,w) = -0{<f>, (3.99) 

The triple summation 9(<p, tp,co) = 0 if (j) = xp or (f> = u, that is, if n — 0 or n — 2, 
respectively. Then, for v = 0, (3.97) indicates that the following (and only the following) 
are ideal invariants for 2-D homogeneous Euler turbulence: 


Ek = l £* 2 |if(k )| 2 = J £ |u(k)|’ 

z kes z kes 

Q = l E^IWI 2 = -El^MI 2 - 

Z keS Z keS 


Again, Ek is the kinetic energy and fi is the enstrophv. 


(3.100) 

(3.101) 


Kelvin’s Theorem 


The x-space equation corresponding to (3.96) with v = 0 is 


doj 

dt 


+ u • Vuj = 


du> 

dt 


0. 


(3.102) 
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This equation tells us that the vorticity of each fluid element is conserved in 2-D ideal flow as 
it is convected about the flow region. This is Kelvin’s circulation theorem. Since duj/dt = 0, 
we also have 


c o 


a— l 


du 

dt 


dco a 

dt 


0 , a>l, 


(3.103) 


so that any non-negative power of ui associated with a fluid element is also conserved as 
that element moves around. (We require a > 1 so that (3.103) is defined for all finite co, 
including co = 0.) Thus, the integral of u a over the flow region (here, a periodic box), is also 
conserved, for any a > 1: 


I a = JJ u a dxdy — const. (3.104) 

If a is drawn from the set of real numbers greater than or equal to one, (3.104) appears 
to imply that there are an uncountable infinity of integral invariants for 2-D Euler flow. 
However, if we restrict ourselves to a periodic domain, then a must be a positive integer, so 
that the number of invariants is, in fact, denumerable. As we have just seen, only the case 
with a = 2 actually leads to a non-trivial invariant in a finite Fourier representation. Such 
a quantity is often referred to as a ‘rugged invariant.’ 

The reason that Kelvin’s theorem does not appear to carry over is three-fold. First, 
Kelvin’s theorem essentially states that the vorticity of each fluid point, and thus all powers of 
vorticity of each fluid point, are conserved in a 2-D ideal fluid. This is a statement concerning 
the continuum and is not commensurate with the fundamental nature of a discrete Fourier 
representation. Second, Parseval’s identity guarantees only the equality of Q, as defined in 
k-space by (3.101), and I 2 , as defined in x-space by (3.104). Third, as will be seen in the next 
chapter, the u>(k), in a numerical simulation, essentially comprise a set of random variables 
with a normal probability distribution. For such a set of random variables, the expectation 
values (hn-i) = 0, while (I 2n ) = (2 n- 1)!! (/ 2 ) n , for n — 1,2,3,... . Also, as n -* oo, 
fluctuations in I 2n are found to grow exponentially. Thus, only Ll = I 2 can be expected to 
behave like an invariant, either analytically or computationally - a truly rugged invariant. 
In a sense, the ( I 2n ), n > 1, are also ‘invariant,’ since they are proportional to ( I 2 ) n , but 
their fluctuations are so great that they cannot be considered ‘rugged.’ 


2 -D Ideal MHD Invariants 

In order to determine the integral invariants for 2-D ideal MHD, multiply (3.44) by </>j(k), 
multiply the complex conjugate of (3.45) by </> 2 (k), add the resulting equations, and sum 
over k e to get 


E 

kes L 


~ doj( k) ~ da*( k) 

«<k) ; ( '+«k) ; 


l 


[0(</>i,^,w) + 0(01,7, a) + 0{(j) 2 ^,a)] 
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+ £ fk • B 0 [j( k) <^j(k) - <&>(k) ^*(k)] 

kes 

— £ k 2 j^u;(k) <^i(k) + r/^2(k) a*(k)j. (3.105) 

keS 

A choice of functions (f>\ and will produce an invariant only if the left side of this equation 
is a perfect differential and if the right side vanishes for v = r\ — 0. 

There are only three choices: 1) 4>i = 4>2 = j] 2) 4>\ — a, 4>2 = w; and 3) (j > i = 0, fa = a. 

The first choice (c/q = ip, <p 2 = j ) yields (3.87), so that the total energy E = Ek + Em is an 
ideal invariant for 2-D MHD. Here, E k and El are defined by (3.100) and (3.101), while E m 
and J are defined by 

Em = ll> 2 |4(k)l 2 = «Elb(k)l 2 (3106) 

z kes z kes 

j = J £fc 4 |a(k)| 2 = l £ li(k)| 2 . (3.107) 

Z kgS z keS 

The second choice (<f>i — a, (j ) 2 = to) yields (3.93), so that the cross helicity He is also an 
ideal invariant for 2-D MHD: 

He = J 5>-<k)«(k) = ' £ ir(k)j(k). (3.108) 

z kes z kes 

The third choice (</q = 0, <f >2 = a) yields 


dA 

dt 


£ik-B 0 a*(k)&(k)-2 iiE m , 

kes 


(3.109) 


where A is the mean square magnetic potential: 

A = I £ |a(k)| 2 . 

z kes 


(3.110) 


It is clear from (3.109) that A is an ideal invariant in 2-D MHD, but only if B 0 = 0. 

To summarize the results for ideal 2-D MHD, E, He and A are integral invariants when 
B 0 = 0, but only E and He are ideal invariants when B 0 ^ 0. We now have at hand the 
integral invariants for the six different types of ideal, incompressible, homogeneous turbu- 
lence, and we will use these in the development of the associated statistical theory. Again, 
it must be stressed that these integral invariants are ‘constant’ only to machine accuracy 
during a numerical simulation. Although fluctuations may be slight, they allow the develop- 
ment of a useful statistical mechanics, which describes the behavior of computer models of 
ideal homogeneous turbulence, and hopefully has some relevance to more realistic models of 
turbulence and ultimately to real-world turbulence itself. It is the basics of this statistical 
mechanics that we develop in the next chapter. 
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3.5 References for Further Reading 

Discrete Fourier transforms are covered in detail in [Hamming 86] and in [Bracewell 86]. 

The energy is a well-known classical invariant for any conservative physical system. Other 
invariants for ideal homogeneous fluid and magnetofluid turbulence were recognized more 
recently: O [Kraichnan 75], Hk [Betchov 61], He [Woltjer 58], H M [Elsasser 56], and A 
[Fyfe 76]. 
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Chapter 4 

Statistical Mechanics 


The evolution of the Fourier modes introduced in the last chapter is determined by sets of 
first-order, nonlinear ordinary differential equations which couple components of different 
wave vectors k to one another. The coupling is highly nonlinear, the various Fourier com- 
ponents are observed to behave stochastically and so must be considered random variables. 
This suggests that a statistical description of the multidimensional dynamics of turbulence 
is required. Although a statistical theory has been elusive for real, dissipative turbulence, it 
is possible for computer models of ideal turbulence and follows rather straightforwardly from 
the principles of classical equilibrium statistical physics. The novelty of ideal turbulent flows 
is that several integral invariants exist and must be utilized, whereas in classical statistics 
only the energy is needed. 

Our particular goal is to develop a statistical mechanics that allows us to analytically 
predict expectation values of moments of random variables and to compare these predictions 
with the results of a numerical solution of the equations of motion. In fact, our statistical 
theory must take into account the approximations inherent in computer simulation, specif- 
ically fluctuations during a numerical solution due to time-step size and round-off error. 
These fluctuations actually bring the statistical behavior of a computer model closer to that 
of a realistic physical system. A mesoscale physical system may be generally thought of 
as a small part of a larger ‘heat bath,’ and a computer model may similarly be viewed as 
a small system embedded in a computer, which serves as its heat bath. In equilibrium, 
the interaction between a ‘small system’ and ‘heat bath’ is only through fluctuations, and 
the statistics of such a small system are best described in terms of canonical , rather than 
microcanonical ensembles. In addition to being the most appropriate, a canonical ensemble 
distribution function can be more readily used for actual calculations. Historically, the use 
of canonical ensembles to give a statistical explanation of ideal homogeneous turbulence was 
given the name absolute equilibrium ensemble theory. It is this theory that we develop in 
this chapter. 
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4.1 Turbulence as a Dynamical System 

Each of the independent real and imaginary parts of the velocity and magnetic field com- 
ponents u(k) and B(k) can be associated with the components of a single, time-dependent 
vector v(f) in a higher dimensional space, called a phase space and denoted by T. The phase 
vector v defines the location of a phase point, which completely describes the instantaneous 
state of the system under investigation. (Although this could have been done using x-space 
grid point values, along with the equations of motion in finite difference form, here we choose 
to work in k-space, as the essential dimension of the phase space is manifest, among other 
benefits.) In the last chapter, we transformed the equations of motion of a magnetofluid, 
from x-space into k-space. This effectively created a nonlinear dynamical system to serve as 
a mathematical model of a magnetofluid. 

For homogeneous turbulence, the dimension n d of the dynamical system follows from 
(3.63) for 2-D flows, and from (3.64) for 3-D flows. In 2-D, we use the complex scalar 
functions d)(k) and a(k) to define the phase space, while in 3-D, we use three complex 
components of each of the vectors u(k) and b(k) to define the phase space. However, in 3- 
D, the three components are kinematically linked by solenoidality (3.12), so that only two are 
independent. Thus, not all of the phase space axes for 3-D turbulence are independent and 
the number of phase space dimensions is larger than n d , but this redundancy is tolerable since 
there is no inherent reason to prefer any two of the components over a third. The dimension of 
the phase space is determined by grid size, for either 2-D or 3-D turbulence and also whether 
B = 0 or not. In the 2-D Euler case, the dimension of the phase space is n = n d = F 2 , while 
in the 2-D MHD case, the dimension is n = n d = 2 F 2 , where F 2 is given by (3.63). In the 
3-D Euler case, the dimension of the phase space is n = 3 over2n d = 3 over2F$, while in 
the 3-D MHD case, the dimension is n — 3 over2n d = 3F$, where Fz is given by (3.64). In 
planning a numerical implementation, the argument can be reversed, since computational 
speed and available core memory dictate the value of the grid size N D , thereby setting the 
largest practical n and n d . 

The equations of the dynamical systems created using isotropically truncated Fourier 
expansions all have the following appearance: 

d dt = ~ KjVj ' t 4 - 1 ) 

(Above, there is no sum over the repeated indices.) The n-dimensional vector v has com- 
ponents Vj, j = l,...,n, each of which represent a different independent component, for 
example, drawn from the 3-D Fourier vectors u(k) and b(k), or from the 2-D scalars u;(k) 
and a(k). The constant Kj is either ok 2 or rjk 2 , as appropriate for the particular Vj. The 
term fj contains the nonlinear coupling of the different Vj. Thus, (4.1) can represent any of 
the 2-D or 3-D evolution equations discussed in the last chapter. Since (4.1) has no external 
driving force on its right side, it is termed autonomous, so that we are dealing here with an 
autonomous, nonlinear dynamical system. 
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4.2 Canonical Ensembles 

Thus far, we have abstracted the various 2-D and 3-D finite Fourier models of fluid and 
magnetofluid turbulence that appeared in the last chapter into the form of an n d -dimensional 
dynamical system contained in an n dimensional phase space and represented by a phase 
vector v that satisfies (4.1). Effectively, we have the whole system identified with a single 
point in the n-dimensional phase space T (n 1), whose trajectory v(f), we assume, follows 
some phase path , beginning at some initial point v(0), and satisfies certain constraints. Here 
we will specify some of these constraints by initially assigning specific, self-consistent values 
to the ideal integral invariants associated with whichever one of the six different cases of 
homogeneous turbulence we are studying. Each ideal invariant is essentially an equation 
that defines a hypersurface of dimension n — 1 in T. In addition, the solenoidal constraints 
(3.12) define a host of other hypersurfaces for cases of 3-D turbulence. The intersection of 
all hypersurfaces is the invariant subspace or manifold for the particular case being studied. 

If the integral invariants were strictly conserved and known to infinite accuracy, the 
invariant manifold would have either na — 2 or — 3 independent dimensions, depending 
on whether the dynamical system has two or three ideal integral invariants. In the case of 
strict conservation, the probability that a phase point was on the invariant manifold would 
be unity and the probability that it was off the invariant manifold would be zero. 

As an example, in 3-D ideal Euler turbulence the invariants are the energy E = Ek 
(3.76) and the kinetic helicity Hk (3.81). If these were exactly conserved, they would always 
be equal to their initial values: E = and H K = H$. Then, the probability distribution 
would be given by 

Dp = 6(E-EM)5(H k -H [ ° ] ). (4.2) 

The phase function Dp is the microcanonical ensemble probability distribution (or density). 
Again, its use would be required if the ideal invariants were strictly conserved, although its 
use in the calculation of expectation values would be a great challenge, since the integration 
must take place on a rather complicated hypersurface in T (called the invariant manifold). 

However, on a computer, invariant values are defined only to a finite level of accuracy 
and, as the system evolves, errors in the computed values of u(k) and b(k) will occur and 
cause the value of each ideal invariant to be satisfied only up to small fluctuations. This, 
in turn, will cause what would have been the microcanonical invariant manifold to move 
about during a simulation. The result is that the fluctuating phase point sweeps out a 
‘fuzzy’ hypersurface in T. In this case, the distribution of phase points is best described by 
a smooth phase function D (but one that is highly peaked on the average position of the 
‘fuzzy’ hypersurface). Such a phase function is called a canonical probability distribution or 
density. 

Since n » 1, the ‘invariant manifold’ is generally nontrivial and there are many choices for 
its initial points v(0). In fact, on a digital computer, there will be a finite, though very large 
number, of possible v(0), which form the set Vo of all initial values of v at t — 0 consistent 
with the initial values of the ideal invariants and the inherent accuracy of the computer. 
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Again, Vo would be the invariant manifold for a microcanonical system, since it contains all 
v € T that satisf)^ all constraints (to machine accuracy). On the other hand, a canonical 
ensemble is the set V of all v(t), as defined by numerical implementation of (4.1), following 
from all of the possible v(0). Although all v(0) G Vo are equally probable initial values, the 
points v(t) 6 V~ T are not all equally probable, but instead have a probability determined 
by the canonical probability density D. Therefore Vo is the initial manifold, which broadens, 
due to fluctuations, into the (fuzzy) invariant manifold for a canonical system. The form of 
D will be determined after a discussion of Liouville’s theorem. 


4.3 Liouville’s Theorem 


The finite set Vo 3 v(0) can be thought of as a ‘gas’ of noninteracting points which move 
about according to (4.1). If the distribution of these points in T at any time t is given 
by D, then the number of points 8n in the small volume of phase space 5r is 8n = D <5r. 
In analogy with our development of the continuity equation (2.7), we require that the <5n 
remains constant ( i.e ., the number of points is conserved): 


dSn 

dt 


dDST 

dt 

dD 

dt 

0. 


5T 


+ D 


dST 

dt 


The small phase space volume 5T can be expressed as 

8T = n 6 Vi, 

i = 1 


so that 


d5T 

dt 


h dt h 

n civ 
3=1 &3 

E 5 [/?( v ) - K i v i\ II 6v 

j =1 i^j 


E 

j=l L 


dvj 


Kd 


Svj JI Svi 

i+3 


- •£ % «r. 

3 = 1 


(4.3) 


(4.4) 


(4.5) 
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The last step follows from (4.4) and because fj(v) does not contain v 3 (as can be seen by 
looking at any of the evolution equations of the last chapter), which implies 


dfj(v) 

dvj 


0 . 


(4.6) 


Defining k = Kj > 0, it is clear that (4.5) leads to 

<5T = <5T^ e~ Kt and D = D [0] e +Kt , (4.7) 


where <5r^ and are the values of <5T and D at t = 0. The relations (4.7) pertain to a 
phase space volume element which moves along with the points it contains, and the increase 
in the density of points D is balanced by the decrease in <5r, so that 6n = D6T remains 
constant. 

If v = rj = 0, then Kj — 0 — * k = 0, and the dynamical system is conservative, so that 
(4.1) represents ideal homogeneous turbulent flow. In this case, equation (4.6) is a critical 
detail for ensuring that both D and 5r remain constant during ideal flow in phase space. The 
result D = constant is sometimes termed a detailed Liouville theorem , to separate it from 
the Liouville theorem associated with conservative Hamiltonian systems, where Hamilton’s 
canonical equations ensure the constancy of probability density and phase space volume. 

It must be emphasized again that ‘constant’ for a canonical ensemble means ‘constant 
to within small fluctuations.’ Thus, the various ideal integral invariants fluctuate slightly 
about their average values and, for example, the expectation value of fluctuations in energy, 
i.e., the variance of energy, is (( E — (E)) 2 ) ^ 0. The expectation value (Q) of any phase 
function Q is defined as 

(Q) = j QDdT. ( 4 . 8 ) 

Here, the integration is over the whole of phase space — oo < Vj < +oo and D is normalized 
so that (1) = 1. Expectation values will be more fully developed in the following sections. 


4.4 Canonical Probability Density 

The canonical probability density function (PDF) D is a phase function, that is, a function of 
the phase variables Vj in the phase space T. The result that D is constant (again, to within 
canonical fluctuations) for ideal flow means that D is itself a function of other constants 
of the motion, i.e., the ideal invariants. Furthermore, for any system, the phase space can 
generally be decomposed into a Cartesian product of subspaces: 

r = r {1) 0 r (2) 

= v (1) 0 v (2) , v e r, v (1) e r (1) , v (2) <e r (2) . 


V 


(4.9) 
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(It is possible that T^ and may be further decomposed, the process continuing until 
some maximal irreducible set of product subspaces is reached.) Each subspace T^, i = 1,2, 
will then have its own statistical distribution function D^\ i — 1,2, such that 

D = D (1 >Z) (2) . (4.10) 

This relation holds whether or not the separate variables in and interact. 

If the r«, i — 1,2, did not interact with one another, the D^ l \ i — 1,2, will each be 
constant functions of ideal invariants pertaining only to their variables, i.e., their inherent 
summations will be restricted to the smaller set of variables found in either or T^ 2 \ 
rather than the whole set found in T. However, if the do interact to form a quasi-closed 
system, this will reach equilibrium and D will be a constant function of the ideal invariants 
associated with the whole phase space T, although the D^, i — 1,2, will individually no 
longer be constants of the motion. 

To illustrate this point, consider 3-D Euler turbulence, where the invariants are energy 
E — Ek (3.76) and kinetic helicity H K (3.81). In this case, we can form the TW, i = 1,2, 
by choosing two disjoint sets of wave vectors k. Let /C be the following set: 

/C — {k | k m in C k — (4.1l) 

and from this choose two subsets K^ l \ i = 1,2, such that 

K = /C (1) U /C (2) and /C (1) n /C (2) = 0. (4.12) 

Usually the choice k m in = 1 is made, but other values can also be chosen. (Note that k m i n = 1 
rather than k min = 0 is used since all modes associated with k = 0 do not dynamically evolve, 
as the results in the previous chapter show. They are always initialized with zero value and 
retain this value during a numerical simulation.) 

As a specific example, since each k has integer coefficients, two subsets that tend to equal 
size in the limit of large k m ax and which satisfy (4.12) are 

^ — }k | kmin — k ^ k max , k is odd} 

/C (2) = {k I k min < k < kmax, k 2 is even}. (4.13) 

Other examples can be devised, with the only requirement being that the Fourier modes 
associated with a given subset may be dynamically coupled to each other, but can be dy- 
namically decoupled from those associated with the other subset. The example embodied in 
(4.13) has some interesting features, as will be discussed in the next chapter. 

The variables in T^, i — 1,2, are the independent components of u(k), k G KM\ i = 1,2. 

For 3-D Euler turbulence, the energies and helicities associated with T and the T^, i — 1,2, 

are 

E = £ (1) +E {2 \ E {i) = 1 £ |u(k)| 2 , * = 1,2 

ke/CM 

H k = 4 !> + Hf, H<$=\ Y. <= 1 , 2 - 

z k 6JC(0 


( 414 ) 
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(All three components of u(k) and u>(k) are included in the sums above, not just the indepen- 
dent ones.) The sets /C and K^\ defined by (4.11), (4.12) and (4.13), contain modes k ~ -k 
which are equivalent because of the reality condition (3.14). To reduce this redundancy, 
define 

K! = {k | k e K, k • x' > 0} 

K' {i) = {k | k e K®, kx' > 0} , i = 1,2 

x' = x- (2A; max ) _1 (y + z). (4.15) 

Thus, /C' and /C'^ have half the members of K and KS l \ respectively. Using these results, 
(4.14) becomes 

E {i) = Y |u(k)| 2 and H$ = Y u*(k)-u>(k), * = 1,2. 

ke/c'(>) keJC'M 

(4.16) 

The factor of \ appears in front of the sums (4.14) because these count twice every term in 
the sums (4.16). 

Since the canonical probability densities are functions of the ideal invariants, the only 
nontrivial functions which satisfy (4.10) when T (1) and T (2) are in equilibrium are 

D = ^ exp {-aE - / 3H K ) 

£><■> = ^ exp [-<*£«> - fiHf ] , * = 1,2. (4.17) 

The normalizing functions ^ and are called partition functions: 

Z = J exp {—aE — 0H K ) dT 

Z (l) = j exp (- aE (i) - 0H$) dT«. (4.18) 

In (4.17) and (4.18) the parameters a and (3 are not phase functions, but instead are un- 
determined multipliers , or as they are often called, inverse temperatures. Terms such as 
qk • UR(k) are not included in (4.18) because these are zero to within small fluctuations. 
Also, the phase volume elements in (4.18) are 

dr = J|dr (i) , dr w = PJ dr(k), dr(k) = du R (k)du/(k). 

i ke/C'M 

(4.19) 

There are analogous expressions to (4.19) for the other five cases of ideal homogeneous tur- 
bulence. (The modal subspaces T(k) for the different cases of ideal turbulence are explicitly 
identified in Table 4.2.) 
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If the r« are isolated from one another, then each will have different inverse temperatures 
and f3^. In this case and will have the form 

D {1) = 2® exp [~^ )E{i) ~ P (i )h k } , * = 1,2. 

Z {l) = Jexp(-a^E {i) dr (i) . (4.20) 

Here, in general, ^ a ^ and ft 1 '* ^ However, if the are dynamically connected, 
then they will both equilibrate to the same ‘temperatures,’ a and /?, returning them to the 
forms (4.17) and (4.18). 

It has been stated that the canonical invariants, such as E and Hx, fluctuate about aver- 
age values and are constant only in a probabilistic sense, i.e ., the fluctuations are relatively 
small. The probability densities appearing in (4.17) (and their analogous forms for the other 
five cases of ideal homogeneous turbulence) provide a way to give this statement precise 
meaning. Canonical probability densities will, in fact, be highly peaked at the average val- 
ues of the integral invariants, but it must be noted that E and H K within the exponentials 
of (4.17) can take essentially unbounded values within T. 

Also, in looking at (4.17), we see a general and very useful feature of the canonical 
probability densities of ideal homogeneous turbulence: They are all essentially Gaussian 
probability distributions, i.e., the arguments of the exponentials can always be put (through 
coordinate transformations, if necessary) into the form of negative definite (because of the 
minus sign) symmetric quadratic forms. This fact makes calculation of the expectation values 
of moments of the variables of T a straightforward procedure, as will be seen presently. 


4.5 Partition Functions 


Z = / exp (-ah ~ PE ~ ih) dT 


Case B 0 

2- D Euler 

3- D Euler 
2-D MHD 0 

2- D MHD 1 

3- D MHD 0 
3-D MHD 1 


h 

h 

h 

E 

Cl 

0 

E 

Hx 

0 

E 

H c 

A 

E 

H c 

0 

E 

He 

H m 

E 

H c 

0 


Table 4.1: General form of partition functions. 


4.5. PARTITION FUNCTIONS 


45 


In the last chapter the evolution equations and ideal invariants of the six basic cases of 
ideal homogeneous turbulence were presented. Here, the six cases, their ideal invariants, and 
the general form of their partition functions are given in Table 4.1. In order to demonstrate 
the evaluation of the general integral of a partition function, consider, once again, the 3-D 
Euler case: 

Z = [ exp (-aE - pH K ) dT = l[ Z(k), (4.21) 

J keK' 

where the modal partition function Z( k) is 

Z( k) = J exp [— a|u(k) | 2 - i(3u* (k) • k x u(k)] dT(k) 

= Jexp[-A(k)]dT(k), (4.22) 

and the modal six-dimensional volume element dr(k) is given by (4.19). Also, recall that 
each variable in (4.22) is evaluated from — oo to +oo. 

Expanding u(k) explicitly into real and imaginary parts and the argument A(k) in (4.22) 
becomes 


A(k) = a|u(k)| 2 + /3u*(k) ■ u>(k) 


= a [|u/i(k)| 2 + |u/(k)| 2 ] + 2 / 0u / (k) • k x u R (k). 


The integral in (4.22) can be evaluated by either the coordinate transformation 

0, 


or the alternative one 


U/(k) = u/(k) -I- k x Ufl(k), 
a 


Ufl(k) = Ufi(k) - ^k x U/(k). 


Upon using the transformation (4.24) on (4.23), the integral in (4.22) becomes 


Z(k) = J exp 


a 2 — fPk 2 

a 


|Utf(k)| 2 — a|u' 7 (k)| 2 


dT'(k). 


(4.23) 


(4.24) 


(4.25) 


(4.26) 


Here, dr'(k) = du R (k) du' 7 (k). [Using (4.25) would merely have switched the subscripts 
R ^ I in the dT'(k).] 

In order to evaluate (4.26), the parameters a and (3 must satisfy the following inequalities: 

a > 0, a 2 > P 2 k 2 max . (4.27) 
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Provided these inequalities are met, using the well known result 

Z(a) = J exp (—ax 2 ') dx = ^ (4.28) 

produces 

Z( k) = 7T 3 (ot 2 — /3 2 /c 2 ) -3/2 . (4.29) 

Note that although we have found Z(k), and thus Z by (4.21), the quantities a and /3 
are still undetermined, though (4.27) limits their values. Also, note that the two coordinate 
transformations (4.24) and (4.25) show that D is essentially a Gaussian distribution function 
on the variables u' fi (k) and u)(k). 

The modal distribution functions for the other cases of ideal homogeneous turbulence 
listed in Table 4.1 can be found in the same way as the 3-D Euler case just considered. The 
results of doing so are given in Table 4.2. In this Table, the six separate cases are listed, as 
well as their associated phase variables and modal distribution functions. Using these results 
in the next section will allow us to find various expectation values, since the Z( k) normalize 
their associated canonical probability densities. 

4.6 PDFs for Ideal Turbulence 

The canonical probability density function D that gives the expected phase space distribution 
of phase points is 

D = n £>( k )> J i = £ * = 1,2,3 

ke/c' ke/c' 

D(k) = z f k) / exp [-<*/, (k) - /J/ 2 (k) - 7 /s(k)] (JT(k). 

(4.30) 

The modal partition functions Z( k) are given in Table 4.2 and the modal terms Ij( k), 
j = 1,2,3, are given in Table 4.3. These will be used in the next section to determine 
expectation values of moments of phase space variables. 

4.7 Expectation Values 

Now, we move on to determining expectation values (4.8) for moments of the components 
of Uft(k), u/(k), bft(k) and b/(k). In what follows, the expectation values apply to those 
components that are kinematically nonzero. To be more explicit, when k x = k x x or k^ = k y y 
or k z = k z z ( i.e ., only one component of k is nonzero), then solenoidality leads to ^(k^) = 
bj(kj) = 0, j = x.y.z for all t. In these cases all expectation values are zero; in all other 
case, they are as given below. 
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Z( k) = f{k,a,p, 7 ) 


Case 


r(k) 

Z( k) 

a, 0 , 7 (<5 2 = a 2 - */? 2 ) 

2-D 

_ 

£bj,/(k) 

TT 

a > 0, p > -a/k 2 max 

Euler 



ot/k 2 + ft 

a < 0, P > -a/k 2 min 

3-D 

- 

iW(k) 

7T 3 

a > 0, |^| <C a/ k max 

Euler 



(a 2 -p 2 k 2 )* 


2-D 

0 

£tR,/(k) 

7T 2 

a > 0 

MHD 



<5 2 + aj/k 2 

7 > 0, 5 2 > -a^/k 2 max 
7 < 0, 5 2 > -aj/kU n 

2-D 

1 

^R,/(k) 

TT 2 

a > 0, 5 2 > 0 

MHD 


c«,/(k) 

S 2 

7 = 0 

3-D 

0 

««,/( k ) 

7T 6 

a > 0, 5 2 > 0 

MHD 


bfl,/(k) 

(5 4 - a 2 7 2 /k 2 )* 

I 7 I < kminS 2 /a 

3-D 

1 

ur,i( k) 

7T 6 

a > 0, <5 2 > 0 

MHD 


V/(k) 

<5 6 

7 = 0 


Table 4.2: Modal partition functions for ideal turbulence. 
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D( k) = Z *(k) exp [— aii(k) - @I 2 { k) - 7/3 (k)] 


Case 

h( k) 

/ 2 (k) 

hi k) 

2-D 

Euler 

k~ 2 \uj R (k)\ 2 
+ k~ 2 |u)/(k)| 2 

|w«(k)| 2 + |o)/(k)| 2 

0 

3-D 

Euler 

\MV\ 2 

+ |u/(k)| 2 

2k • tifl(k) x u/(k) 

0 

2-D 

MHD f 

k~ 2 |^(k)| 2 
+ k- 2 \uj I (k)\ 2 

+ k 2 \d R (k )\ 2 

+ k 2 \dr(k )\ 2 

ur( k) dfl(k) 
+ w/(k) o/(k) 

|di?(k)| 2 + |d/(k)| 2 

3-D 

MHD t 

|u fl (k)| 2 
+ |u/(k)| 2 
+ |b/?(k)| 2 
+ |b/(k)| 2 

Ufi(k) • bfi(k) 

+ u/(k) • b/(k) 

2k • a fi (k) x a/(k) 


^When B 0 = 1, 7 = 0. 


Table 4.3: Modal density functions for ideal turbulence. 


First, we note that the following integrals of an odd integrand between symmetric limits 
are zero: 

» OO 


/ OO / v 

rr 2n_1 exp (— ax 2 ) dx = 0, n = 1,2,3,.... 

-OO ' 


(4.31) 


Similar integrals of even powers of x are nonzero, however, and can be found from (4.28): 

/ CO * . 

x 2n exp ( — ax 2 ) dx 

-00 ' ' 

1(a) 

1 

1(a ) , n = 1, 2, 3, — 


= - 
V da , 

(2n — 1)!!. 

(2a)' 

In terms of (4.28) and (4.32), the expectation value of an even power of a variable with a 


(4.32) 
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Gaussian distribution function is 


/„2n\ - ^ 2n ( a ) 

' / 1(a) 

(2n-l)\\ „ 19 , 

(2a) n ’ 1>2,3,.... 

(4.33) 

In particular, for n = 1 , (4.33) gives 



<* 2 > 

II 

H ^ 
II 

S? - 

(4.34) 

while for n > 1 , we have 



V") 

= (2n - 1)!! (x 1 ) 71 . 

(4.35) 


Thus, knowledge of ( x 2 ) gives knowledge of all the nonzero moments (x 2n ), thereby com- 
pletely determining the statistics of x, which is a well known result for Gaussian distributions. 

The expectation values of powers of the phase variables can now be calculated for the six 
cases of ideal homogeneous turbulence using the general form 

(u m ( k)) = Ju m (k)DdT 

= J u m (k)D(k)dT{k), (4.36) 


where the modal canonical probability densities D(k) are defined in (4.30) and given explic- 
itly in Table 4.3. Evaluation of expectation values for the different cases of ideal homogeneous 
turbulence requires the use of coordinate transformations such as (4.24) and (4.25). Since 
these expectation values are important, we will look at them in detail. The general form of 
the probability densities is 


D( k) 


1 

Z(k) 


exp [-A(k)] 


A(k) = cr/j^k) + /3/ 2 (k) + 7 / 3 ^). 


(4.37) 


The values of /;(k), i = 1 , 2 , 3, for the different cases are given in Table 4.3. 


4.7.1 2 -D Euler 

In the 2-D Euler case, the modal argument in (4.37) is summation over the real and imaginary 
parts of a modal variable: 

-4(k) = Y [«^I(k)A 2 +/3D|(k)] 

S=R,I 

= Y (aA 2 + P) w|(k). 

S=R,I 


(4.38) 
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If we put the following into (4.31) and (4.34), 

£ = a>s(k), a = a/k 2 + (3, (4.39) 


we easily determine 

fe(k)) = 0, ("!«)= 2(^*2). S = R,I. (4.40) 

Expectation values of higher even powers can be generated using (4.33). 


4.7.2 3-D Euler 

In the 3-D Euler case, the modal argument in (4.37) is given by (4.23), after applying the 
transformation (4.24): 

A(k) = a [|u H (k )| 2 + |u/(k)| 2 ] + 2f3u I (k) • k x u fl (k) 

= ex' 1 ( a 2 - /3 2 /c 2 ) |u«(k )| 2 + a^u^k)! 2 . (4.41) 

As was noted earlier, (4.25) can be used instead of (4.24), which merely switches R and I in 
(4.41). Thus, using these results, along with (4.31), (4.34) and (4.37) gives 

( , ^s , ,i(k)) 0, S R , 7, i x, y, z , 

«i(k)) ^ 2 (a 2 - p 2 k 2 ) ' ^ 4 ' 42 ^ 


Since (4.41) also gives 

(«j,i( k )ttj,i( k )) = 2a’ i = x ’y> z ’ ( 4 - 43 ) 

we can use (4.24), along with our ability to interchange R and /, to show 

—/3k 2 

(ws,i(k) u>s,i(k)) = 2(cP-^k 2 y S = R ’ R i = x ^V^ z - 

(4.44) 


4.7.3 2-D MHD 

In the 2-D MHD case, the modal argument in (4.37) is 

A(k) = Y, {» [*r 2 D|(k) + & 2 a|(k)] + £w s (k)as(k) + 7a|(k)} 

S=R,I 
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= E 


S=R 


a 


j Ik 2 


U)\ 


;(k) + /8ws(k)a s (k) + ( ak 2 + 7) a|(k) 


= E 

S=R,I L 

= E 

S=R,I L 


"ws 2 ( k ) + * 2 + 7] «K k ) 


^ ) ^(k) + («fc 2 + 7) «5 2 (k) 


(4.45) 


ak 2 + 7 

Here we have introduced 8 2 = a 2 — \(3 2 for brevity. In (4.45), the two transformed variables 


are 


Bk 2 

c^(k) = 65 s (k) + P 2a a 5 (k), S = R,I 
S' s (k) = ^(k)+ 2(a ^ +7) ^(k). 


(4.46) 


From these results, it is a straightforward matter to deduce the following expectation values: 

(£>s(k)> = (os( k)) = 0, S = R, I 
ak 2 + 7 


( w s( k )) 2 (S 2 + a'y/k 2 ) 

(«s( k )) = 

(w s (k)a s (k)) = - 


2 k 2 ( 5 2 + aj/k 2 ) 

P 

4 (8 2 4- a'y/k 2 ) 


(4.47) 


Again, S 2 — a 2 — \P 2 . In the case where B 0 = 0, we have 7 7^ 0, while in the case where 
B 0 = 1, we have 7 = 0. 

4.7.4 3-D MHD 

Finally, in the 3-D MHD case, the modal argument in (4.37) is 

A(k) = Y [ Q! |us( k )| 2 + a|bs( k )| 2 + p\is(k) ■ bs(k)] 


S=R,I 


7 


+2 fc ;k.b#)xb / (k). 


(4.48) 


There are essentially two ways to transform variables in this argument. The first starts with 

Us( k ) = u 5 (k) + ^b 5 (k), S — R,I, 


(4.49) 
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and then uses either of the following 

b',(k) = b,(k)+ ”^kxb,(k) (4.50) 

b' H (k) = b I (k)-^ ! kxb,(k). (4,51) 

Using (4.49) and (4.50) transforms (4.48) into 

x4(k) = alu^^p + alu'^k)! 2 


+ 


5 2 — a 2 j 2 /k 2 


ib R (k)i 2 + ^ibUkjp. 


aS 2 ' ' ' a 

If (4.51) had been used, R and I would be switched in (4.52) 
The second approach uses 

u'/(k) = u/(k) + 2 ^b fi (k), 


followed by (4.50), and then 

b'fl(k) = b*(k) + 


ns 2 


bfi( k )- 


2a(5 2 — j 2 /k 2 ) 

Using these transforms (recalling <5 2 = a 2 — /3 2 /4) turns (4.48) into 

^( k ) = a(^- 7 V / k 2) l “ E(k)|2 + “ l “' <k>|2 


S 2 


a 


(4.52) 


(4.53) 


(4.54) 


(4.55) 


A similar set of transformations would produce an equation identical to (4.55), but with R 
and / switched. 

The expectation values for the 3-D MHD case are thus: 

(u s (k)> = (b s (k)) = 0, S = R,I 


(u s( k )) 
(|bs( k )| 2 ) 

(u s (k) • b s ( k )) 

(a 5 (k)-b 5 (k)) 


3 a(6 2 — 7 2 /A: 2 ) 

2 <5 4 — a 2 7 2 /A; 2 

3 a6 2 

2 5 4 — a 2r y 2 /k 2 

_3 /35 2 

4 S 4 — a 2, y 2 /k 2 

3 a 2r y/k 2 

2 d 4 — a 2 j 2 /k 2 


(4.56) 
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Once again, S 2 = a 2 — */? 2 . In the case where B 0 = 0, we have 7^0, and in the case where 
B 0 = 1, we have 7 = 0. The factor of 3 appears in these results because u(k) and b(k) are 
3- vectors, whose three components all add equally to (4.56). 

4.8 References for Further Reading 

One good reference on dynamical systems is [Verhulst 96]. 

Microcanonical and canonical emsembles are discussed in [Khinchin 49] and [Landau 80]. 

That a detailed Liouville theorem exists for ideal turbulence was discovered by [Lee 52], 
Initial applications of canonical ensemble theory to the six cases of ideal turbulence were: 2- 
D Euler [Kraichnan 75]; 3-D Euler [Moffatt 69, Kraichnan 73]; 2-D MHD ( B 0 = 0) [Fyfe 76], 
(B 0 = 1) [Shebalin 83]; and 3-D MHD ( B 0 = 0) [Frisch 75], (B 0 = 1) [Shebalin 89, Shebalin 94] 
A good review article on 2-D turbulence, including some discussion on ideal 2-D Euler and 
MHD turbulence ( B 0 = 0) is [Kraichnan 80]. 
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Chapter 5 

Thermodynamics 


Canonical ensemble theory, as it applies to ideal homogeneous turbulence, was discussed in 
the previous chapter and allowed us to determine the statistical mechanics of an interacting 
system of Fourier modes. In fact, the independent Fourier modes, denoted by unique wave 
vectors k e K' (see eq. 4.15), form the set of ‘molecules’ of the interacting system, in 
analogy with the individually numbered real molecules of a Maxwellian gas. A gas molecule 
has its own set of associated variables; e.g., in a polyatomic gas, each molecule has its 
own momentum and angular momentum, each vector having three components, for a total 
six variables (neglecting vibration). This defines a 6-D subspace in the phase space of the 
Maxwellian polyatomic gas. The modal subspaces of ideal turbulence are defined in exactly 
the same way, so that the modal subspaces F(k) in the different cases of ideal turbulence 
have their own characteristic dimension (though not all dimensions are independent): 

m = dimT(k). (5.1) 

For the different cases, the modal variables are listed in Table 4.2 and give, for the different 
cases of ideal turbulence, m — 2 for 2-D Euler, m = 4 for 2-D MHD, m = 6 for 3-D Euler, 
and m = 12 for 3-D MHD. Thus, each Fourier mode is a ‘molecule’ with m degrees of 
freedom. 

Furthermore, knowledge of the invariant quantities for the different cases, as summarized 
in Table 4.1, enabled us to find the corresponding canonical probability densities and use 
these to determine expectation values for the moments of the modal subspace variables. This 
is the essential product of the statistical approach - ‘microscopic’ quantities, i.e., modal ex- 
pectation values. What we wish to do now is move on to look at ‘macroscopic’ quantities, 
such as the inverse temperatures a, (3 and 7 (which are heretofore undetermined). These 
macroscopic quantities pertain to all modes and are, in fact, found by summing over all 
independent modes k. This summation over all modes moves us from the realm of statis- 
tical mechanics into the realm of thermodynamics , where we are no longer concerned with 
specific modes, but instead with summations over all modes. These summations produce 
thermodynamic functions that will have specific values for the closed system of interacting 
Fourier modes. 
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Thus, in analogy with the theory of a molecular gas, we move from a statistical description 
to a thermodynamic one. To maintain the analogy, we must ensure that analogous ‘laws of 
thermodynamics’ are present in our model system. The first law is conservation of energy, 
and this is always a invariant for ideal turbulence. The second law, i.e., that entropy never 
decreases after an interaction, will be established in this chapter for ideal turbulence. (The 
third law, the vanishing of entropy as temperature goes to zero, is quantum mechanical, and 
not applicable here.) To begin the process, we first discuss temperature (or rather inverse 
temperatures) . 


5.1 Inverse Temperatures 

We have already encountered several thermodynamic functions, namely, the energy, enstro- 
phy, mean square vector potential, and the various helicities. We can introduce others, as 
required, to facilitate our analysis. All of these thermodynamic functions will have expecta- 
tion values for the system under consideration, which can be determined in the established 
manner. Some of these functions, in particular the invariants already identified, will have 
essentially constant values as the associated system evolves with time, while others, known 
to vary with time, will nonetheless also have specific expectation values. The difference 
between the invariants and the non-invariants is that we can predict the expectation values 
of the first group, but cannot predict those of the second. A priori , then, we must treat 
the expectation values of non-invariants as unknown parameters, a role which will turn out, 
in fact, to be useful. Let us give these preliminary remarks some substance by explicitly 
considering the different case of ideal homogeneous turbulence. 


5.1.1 2-D Euler 


In this case, the two functions we have already met are the energy E and the enstrophy Q: 

E = Y1 ^ _2 |^(k)| 2 (5-2) 

keio 

Q = 22 Mk)| 2 ' (5-3) 

ke?C' 

If we take the expectation values of the above, and use the modal results given in (4.40), we 
get 


(E) = £ 
(f2) = Q 


E 1 

k fe a ^ 

k 2 

hi o 2 ' 


(5.4) 

(5.5) 


Thus, we see that the inverse temperatures a and /3 are implicit functions of £ and Q. While 
these two equations can be solved numerically, the procedure is somewhat problematic (and 
becomes more so when there are three inverse temperatures). 
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Instead, we will look at the expectation value of the non-invariant mean square stream 
function T: 


which is 


$ = k 4 l^(k)| 2 , 

kefc 


^ ^ k 2 (a + pk 2 )' 


(5.6) 


(5.7) 


Although £ and O will differ from the initial values of E and Q, only up to small fluctuations, 
we do not know what 9 is, a priori However, we can use 9 to parameterize a and ft by the 
following procedure. 

First, we note that the following algebraic relations can be derived from (5.4), (5.5) and 
(5-7): 


a£ + PO 

= £ i = v. 

ke/c' 

(5.8) 

a9 + p£ 

= T 1 = c 

h.2 

k 6/C' ^ 

(5.9) 


These ‘thermodynamic relationships’ are simultaneous linear equations and are readily solved 
to yield: 


co - a re 

09-s 2 

Af'9 - C£ 
09 - S' 2 ' 


(5.10) 

(5.11) 


These have been written so that the denominators on the right sides are positive: 


09 -£ 2 


1 (q 2 - k 2 ) 2 

2 qMeic ( a + ^ 2 )( a + &k 2 ) 


> 0. 


The numerator in (5.10) is: 

co - a re 


a 

2 


£ 

q,kE/C' 


(q 2 — k 2 ) 2 

q 2 k 2 (a + Pq 2 )(a + fik 2 ) 


The numerator in (5.11) is: 

A f'9 - C£ 


p ( q 2 _ fc 2 y 


(5.12) 


(5.13) 


(5.14) 
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Placing (5.12), (5.13) and (5.14) into (5.10) and (5.11) yields a tautology. 

At this point, we can allow & to vary: 'T — > which turns (5.10) and (5.11) into 




CQ - APS 
fW - £ 2 

A rw - C£ 
w-s 2 ' 


(5.15) 

(5.16) 


Here, £ and fl are constant values and W is a variable parameter, > £ 2 /f2, such that 
ol = a and /?' = /?, when = ty. The challenge is to find a procedure for determining when 
T' — . so that a and (3 may be found, and the modal expectation values (4.40) exactly 

determined. This will be done in the next section, after we treat the other cases of ideal 
turbulence and develop results analogous to the above expressions. 


5.1.2 3-D Euler 

In this case, the two (invariant) functions are the energy E and the kinetic helicity H K : 


E = T. l“( k )P 

kex:' 

(5.17) 

Hr = £u(k)-w(k). 

k€/C' 

(5.18) 

If we take the expectation values of these, and use the modal results 
(4.44), we get 

given in (4.42) and 

{E) = £ = k g, 

(5.19) 

(H K ) = Hk = Y* a 2_02 k 2- 

(5.20) 

We see that in the 3-D Euler case, the inverse temperatures a and /3 are implicit functions 
of £ and Hr- 

Here, we will look at the expectation value of the enstrophy 0, which is non-invariant for 
3-D Euler turbulence. The enstrophy is 

D = ^2 k 2 |u(k)| 2 , 

keic 

(5.21) 

and its expectation value, using (4.42), is 


(9) = n = y 3k2a 

(5.22) 
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Both £ and Hk will differ little from the initial values of E and Hk , but in this case we do 
not know the value of T2, a priori. Here, we will use a variable f2' to parameterize a and /?, 
in a similar manner as before. 

Using AT as defined in (5.8), the following algebraic relations can easily be derived from 
(5.19), (5.20), and (5.22): 


a£ + £Hk = 3 AT 
aU K + pn = 0. 


(5.23) 

(5.24) 

These simple linear equations readily provide a and /3: 



3API2 

a ~ £Q-H\ 


(5.25) 

n Hk 

13 = - ft “• 


(5.26) 

Again, we have a thermodynamic relationship: the inverse temperatures for 3-D Euler tur- 
bulence are given explicitly in terms of the unknown value Q. Thus, the equations 

, 3 Af'O' 

a ~ ££2' - Hi 


(5.27) 

"a 

1 

II 


(5.28) 

produce the correct values a and ft when the variable parameter f2' is set equal to 12, 
in turn allows the modal expectations values (4.42) and (4.44) to be determined. 

which 


5.1.3 2-D MHD 


Here, the (invariant) functions are the total energy E, the cross helicity He, and the mean 
square vector potential A (for B 0 = 0): 

E = Yl [& -2 |o>(k) | 2 + k 2 |a(k)| 2 ] (5.29) 

ke/c' 

H c = £<h*(k)a(k) (5.30) 

ke/c' 

a = Y K k )l 2 - ( 5 - 31 ) 

ke/c 1 


If we take the expectation values of these, and use the modal results given in (4.47), we 
arrive at 


(E) = £ 


£ 

k eK' 


2a + 7 /k 2 
5 2 + 0:7 / k ? 


(5.32) 
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(He) 

(A) 


He = Y ^ 2 
5 2 + c * 7 /* 2 

A = Y a ^ 2 

lew <5 2 + «7A 2 ' 


(5.33) 

(5.34) 


Here, 8 2 = a 2 — ]/3 2 . We see that in the 2-D MHD case, the inverse temperatures a. f3, and 
7 are implicit functions of (E), (He) and (A). 

In this case, we also consider the expectation value of the magnetic energy Em, which is 
a non-invariant function for 2-D MHD turbulence: 


E m = £fc 2 |d(k)| 2 

ke/C' 

The expectation value of this, using (4.47), is 

(Em) = £m = y 


a 


(5.35) 


(5.36) 


keJC , 5 2 + aj/k 2 ' 

Again, 5 2 — a 2 — \ j3 2 . While £. He, and A are constants, the value of £m will have relatively 
large fluctuations. Here, we will use £' M to parameterize cc, /?, and 7. 

Using AT as defined in (5.8), the following algebraic relations can easily be derived from 
the above expressions: 


a£ 

II 

+ 

0 

£ 

+ 

2 AT 

(5.37) 


2 aHc + P£m — 

0 

(5.38) 

a(£ 

— 2 £m) — J-A. — 

0. 

(5.39) 

s’ readily provide a, /3, 

and 7: 


a 

N'£ m 

£m(S — Em) — 

H 2 C 

(5.40) 

P 

s 

0 s 

£ <0 
1 

II 


(5.41) 

7 

£ — 2£m 

" A “■ 


(5.42) 


The inverse temperatures for 2-D MHD turbulence are thus given explicitly in terms of the 
unknown value of Hm- Here, we use a variable £' M to parameterize the inverse temperatures: 

X'S'm 

£'m(£ - £' u ) - -Hi 


a = 


nr / 

P = 

C M 


(5.43) 

(5.44) 


7 = 


£ - 2£' a 

A 


M a 1 . 


(5.45) 
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Again, the true value £' M = £m is unknown, but once it is found, allows us to determine 
a, /?, and 7 and thus the exact values of the modal expectation values (4.47). (In the case 
where B 0 = 1, we require 7 = 7 ' = 0, which leads to £' M = \ £ as the correct value, without 
further ado.) 


5.1.4 3-D MHD 


Here, the (invariant) functions are the total energy E, the cross 
square vector potential H M (for B 0 = 0 ): 

helicity He and the mean 

E = £ [|u(k)| 2 +|b(k)| 2 ' 

ke/c # 

(5.46) 

H c = £u‘(k)-b(k) 

ke/c' 

(5.47) 

E m = £ a*(k)-b(k). 

ke/c' 

(5.48) 

If we take the expectation values of these, and use the modal results given in (4.47), we 
arrive at 

m c V- 3o;(2(5 2 - 7 2 //c 2 ) 

( } “ 5 4 — a 2r ) 2 /k 2 

(5.49) 

(tt \ _ n, _ Y' 3/3A 2 /2 

' Hc ' ^.,S 4 -a 2 ^/k 2 

(5.50) 

(h \ - u - v 

(^) Em ]L' S 4_ a 2 y 2/k2- 

(5.51) 

Here, 5 2 — 0 ? — \£ 2 . We see that in the 3-D MHD case, the inverse temperatures a, £ and 
7 are implicit functions of £, Ec, and Em- 

As in the 2 -D MHD case, we utilize the expectation value of the magnetic energy E M , 
which is also a non-invariant function for 3-D MHD turbulence: 

Em = E^IMk)| 2 . 

k eK 1 

(5.52) 

The expectation value of this, using (4.47), is 


(Em) — £m = ^2 M 2 2/t2' 

h 4 ~ a?Y/k 2 

(5.53) 

Again, 5 2 = a 2 — \/3 2 . While £, Ec, and Em are constants, 
relatively large fluctuations. Here, as in the 2-D case, we use £' M 

the value of £m will have 
to parameterize a, ft, and 


7 - 


62 


CHAPTER 5. THERMODYNAMICS 


Using A P as defined in (5.8), the following algebraic relations follow from the above 
expressions: 


olE + jSHc + 

- 6AT 

(5.54) 

2 orHc + /3£m 

= 0 

(5.55) 

a{£ — 2 E m ) — i'Hm 

= 0. 

(5.56) 


These ‘thermodynamic relations’ (which are the same as the 2-D MHD case if Af and A 
there are replaced by 3A f and Hu here) immediately provide a, /?, and 7: 


ZN’Em 

“ = £ u (e-s M )-n 2 c 

(5.57) 

0 = 

(5.58) 

£ — 2£ m 

T = v a. 

tiM 

(5.59) 

The inverse temperatures for 3-D MHD turbulence are given explicitly in 
known value of Hu- Using the variable £’ M to parameterize these gives: 

terms of the un- 

3 M’£' m 

£'m(£ — e' M ) — He 

(5.60) 

-3 

0 

£ <0 

CM 

1 

II 

QCi. 

(5.61) 

y = s 2£ ' M a'. 
Hm 

(5.62) 


Although the true value S' M = £m is unknown, once it is found, allows us to determine a, (3, 
and 7 and thus the exact values of the modal expectation values (4.56). (In the case where 
B 0 = 1, we require 7 = 7' = 0, which again leads to £' M = \ £ as the correct value.) 


5.2 Entropy 

Now that we have expressed the inverse temperatures for all cases of ideal homogeneous 
turbulence in terms of a variable parameter, we can move on to defining a procedure for 
actually determining the expectation value of the parameter. We do this by first defining 
the entropy for each of the cases of ideal turbulence, and by showing that the entropy we 
define in each case satisfies the second law of thermodynamics. The probability density for 
all cases has the form 

D = Z - 1 exp (-aR - /3I 2 - 7/3) (5.63) 

z = n 2 < k ). 

kefC' 


(5.64) 
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where the Ij , j = 1,2,3, are given in Table 4.1 and the Z( k) are given in Table 4.2. The 
entropy is defined as 

5 = — (In D) 

— ot ( I\ ) + f3 (h) + 7 (h) + InZ 

= cN'+Y, lnZ(k). (5.65) 

ketc 1 

Here, c — 1 m is the ‘specific heat’ of a mode, where m is defined by (5.1) as the number 
of degrees of freedom per mode. The term cJ\P in (5.65) arises when the expressions (5.8), 
(5.23), (5.37) and (5.54) are used for the different cases of ideal turbulence. Thus, c = 1, 2, 3, 6 
for the 2-D Euler, 2-D MHD, 3-D Euler, and 3-D MHD cases, respectively. 

The Z( k) are functions of the inverse temperatures a, f3, 7, and k as Table 4.2 shows. 
Since a, /3, and 7 are implicit functions of the invariants, the entropy 5, as defined in (5.65), 
is also an implicit function of the invariants and thus is also a constant for a closed system; 
i.e., it does not change unless the system is allowed to interact with another previously 
isolated system. As (5.65) shows, entropy is also a function of the number of independent 
modes M' in the set of modes /C under consideration. (The subset of independent k is 1C.) 
The entropy has the following functionality: 

5 = a()C,a,/3, 7). (5.66) 

Using (5.65) as well as Table 4.2 allows us to list the entropy functionals o{K, a, (3, 7) in 
Table 5.1. 

Now, if we use the variable parameters J?', or £' M to define variable inverse tempera- 
tures a', /?', and 7', as we do in the previous section, then we can allow the entropy to vary 
from its expectation value of S to S', where 

S' = a(fC,a',/3',C). ( 5 . 67 ) 

We will now show that a()C,a',j3', 7') takes a minimum value when a' = a, ft — ft, and 
y = 7. This will then provide us with a procedure for finding a , /3, and 7 without having to 
solve multivariable implicit equations, and more importantly, it will enable a demonstration 
that entropy never decreases when two previously isolated systems are brought together into 
one interacting system. 

5.2.1 A Global Minimum 

As a specific example of the global behaviour of entropy with respect to the variable param- 
eters that have been defined, we will examine the 2-D Euler case. Using Table 5.1, equation 
(5.67) becomes 

S' = A/" r (l -h In 7r) — \n(a'/k 2 + /3') 

keic' 

— N'{1 + ln7r) + y. In A; 2 — y* t In (a 1 + (3'k 2 ). 
keic k eK' 


(5.68) 
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a()C, a, p,j) = c[Af'(l + In 7 r) - Y.k> ln «, P, 7)] 


Case 

c 

(f>(k,a,p,y) 

2-D Euler 

1 

a/k 2 + P 

3-D Euler 

3 

(a 2 — p 2 k 2 ) 2 

2-D MHD f 

2 

(5 2 + a.^/k 2 ) 2 

3-D MHD f 

6 

(5 4 — a 2 7 2 //c 2 ) 4 

tHere, 5 2 = 

a 2 - IP 2 ', 

also, when B 0 = 1, 7 = 0 


Table 5.1: Entropy functionals for ideal turbulence. 

Using a'pE'), as given by (5.15), and P'i'T'), as given by (5.16), we have 

da' _ fta' 
dSP = ~ft<T '-£ 2 

dp' _ Af 1 - Tip' 
d#' ~ ft$'-£ 2 ' 

Using these results, the derivative of S' by is 

dS' _ OS' da' dS'dP' 
dW ~ da' dW + dp' dtif 

^ 1 a' ftp (p'ft - A f')k 2 

k ~^, a' + P'k 2 ft'T' — £ 2 

_Af'(ft'-ft) 
ft $’-£ 2 ' 

(Recall that \P' > £ 2 /ft.) Here we have defined ft' as 

k 2 

ft' — y' 

k «' Q '+^'* 2 


(5.69) 


(5.70) 


(5.71) 
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The result (5.70) indicates that dS' /d,F' — 0 when FI' — FI = 0, or 

k 2 k 2 


fl'-FI = Y. 

k eK' L 


£ 

ke^c 


a' + /3'k 2 a + /3k 2 

(a' - a)k 2 + (£' - £)A; 4 
(ex' + /3'k 2 )(a + /3k 2 ) 


(5.72) 


Now, we can use (5.10), (5.11), (5.15), and (5.16), along with some algebraic manipulation, 
to show 


ex' — a = -(/3'-/3)k 2 , k 2 = 

Placing this into (5.72) gives 

FT -FI = (p'-/3)F' 

k 2 (k 2 -k 2 ) 

Gt K , (a' + /3'k 2 ) (a + /3 k 2 )' 


(5.73) 


(5.74) 


We know from (5.70) that S' has an extremum at FI' = FI, and we now prove two assertions: 
1) that this only occurs at F' = F, and 2) that the extremum is actually a minimum. 

To prove the first assertion, we must show that F' ^ 0 when Fl' — Fl (initially leaving 
open the possibility that F' ^ F). Using (5.10), (5.11), (5.15), and (5.16) again, along with 
some more algebraic manipulation, shows that 


F = F'\n>=n = 


a'£ 


(k 2 - k 2 f 


Af' 2 k ^, (a' + f3'k 2 )(a + /3k 2 ) 

Here we obviously have F > 0, so that F^0. Using (5.11) and (5.16) gives 


(5.75) 




£a(F' - F) 
FIF'-£ 2 ' 


(5.76) 


Placing (5.76) into (5.74), as FI' — > FI, we also have F' — > F, so that 


lim FI' = FI-(F' - F)G , 


where G , using (5.75), is defined by 

G = 

Using (5.70) and (5.77), we see that 


£ 2 a ^ ( k 2 - k 2 ) 2 

A f' 2 (FIF - £ 2 ) (a + /3k 2 ) 2 


(5.77) 


(5.78) 


,. dS' A f'G lTl T . 


(5.79) 
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Thus, we have proven our first assertion, that an extremum occurs only at 'T' = T. 
Next, using (5.79), we immediately see that 


d 2 S' _ A f'G 

d*' 2 ” W-S 2 


(5.80) 


This proves our second assertion: that = & is a minimum of S', and by the results given 
above, it is a global minimum. 

Although it will not be shown here, we assume that analogous results for the other cases of 
ideal homogeneous turbulence can be found, requiring only increasingly laborious algebraic 
manipulations. (These manipulations have been carried through for the 3-D Euler case, 
although not for the ideal MHD cases with B 0 = 0. It is conjectured, and fully expected, 
that the MHD cases will follow suit.) The results, in each case, should mirror what has been 
shown above, that when the entropies as defined in Table 5.1, and the inverse temperatures 
are given a variable parameter, as described earlier in this chapter, then the entropies have 
a global minimum with respect to this parameter. Now, we will move from these facts and 
plausible conjectures to establish a ‘second law of thermodynamics’ for ideal homogeneous 
turbulence. 


5.3 The Second Law 

First, a heuristic look at the expansion of a restricted set of modes into a larger set will be 
given, and second, a general proof that entropy never decreases when two previously isolated 
subsystems are brought together - The Second Law of Thermodynamics. 

5.3.1 Expansion Into a Larger Set of Modes 

To begin, we define two noninteracting systems that can then be allowed to interact. We 
use the two disjoint sets of modes KS l \ i = 1,2, defined in (4.13). Recall that both these 
sets consist of vectors k, 0 < k < k max , with the difference between them being that Kf 1 ^ 
has only those k with k 2 odd, while K. ^ has only those k with k 2 even. 

Now, k 2 is the sum of two squared integers (‘squares’) for 2-D flows, while it is the sum 
of three squares for 3-D flows. The following results from number theory are applicable here: 

Sums of two squares: k 2 is the sum of two squares if and only if its (unique) prime factor- 
ization contains no odd powers of primes p, such that p = 4m + 3, where m is a nonnegative 
integer. 

Sums of three squares: k 2 is the sum of three squares only if it is not of the form 4” (8m + 7), 
where n and m are nonnegative integers. 

Thus, k 2 cannot take all possible positive integer values, in either 2-D or 3-D. (There is 
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a well-known theorem, however, that states that sums of four squares can take all nonnega- 
tive integer values.) 

The number of positive integer values that k 2 < can take out of the total number 
available ( i.e ., k 2 nax ) is different for 2-D and 3-D. A simple computer program shows that, in 
the 2-D case, most positive integers < A;^ aa . are not the sum of two squares (for = 10 4 , 
k 2 takes 2,749 out of 10 4 possible values - about 27.5%; for = 10 6 , k 2 takes 216,314 
out of 10 6 possible values - about 21.6%). In the 3-D case, the ratio of integers that are the 
sum of three squares for a given k 2 max to those that are not is almost exactly 5 to 1 (about 
83%), as k% iax increases. However, in both the 2-D and 3-D cases, the number of individual 
vectors k in the sets i = 1,2, appears to be essentially equal, for k ma x >10. 

The two sets and JC^ have an interesting structure: KS 1 '* identifies a completely 
noninteracting set of wave vectors, while comprises a set of completely interacting wave 
vectors. To see this more clearly, we will now examine the parity (even-oddness) of the 
components of the vectors in /C (1 ) and Let us represent even integers by ‘e’ and odd 

integers by ‘o.’ Adding or multiplying even and even, even and odd, or odd and odd integers 
yields even or odd integers in the following pattern: 

e 4“ 6 ^ 6 0 -|“ O rsj ' 6 6 ~h O = O ~h 6 ^ 0 

(5.81) 

e - e ~ e o • o ~ o e • o — o • e ~ e 

Now, remember that for k G KS l \ we have k 2 ~ o, while for k € KS 2 \ we have k 2 ~ e. Using 
the equivalences in (5.81), it is easy to see that, in 2-D, the kinds of vectors found in the two 
subsets and are 

/C(!) /cW 

(e,o) (e, e) . (5.82) 

(o, e) (o, o) 

Again using (5.81), in 3-D, the kinds of vectors found in K . ^ and 1C ^ are 

/C(!) /C( 2 ) 

(o, o. o) (e, e, e) 

(o,e,e) (e, o, o) . (5.83) 

(e,o,e) (o,e,o) 

(e,e,o) (o,o,e) 

Since a component can have one of two values (e or o), the total number of types of vector 
is N = 2 n , where n is the dimension of the vector space. Thus, for 2-D, N = 4, and all 
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possible types are accounted for in (5.82). Similarly, for 3-D, N = 8, and all possible types 
are accounted for in (5.83). 

We can use these results to deduce the following: When two vectors p, q G are 
added, then p + q = k ^ KP\ because p + q ~ e; thus no p,q 6 can be added to 
produce another wave vector in On the other hand, for p, q G KP\ then p + q = k ~ e; 
if k < k max , then k G KP\ Thus, two p.qe K 'PI can be added together to produce another 
k G KP\ In (3.59) and (3.60) it was shown that any nonzero k G /C could be written as 
the sum of two nonzero wave vectors p, q G /C; looking at the details there, we see that we 
can also state that any nonzero k G KT 2 ' 1 could be written as the sum of two nonzero wave 
vectors p, q G KP\ for > 2. Thus, none of the k G JC ^ are dynamically coupled, while 
all of the k G are coupled. 

The relevance of all this for the simulation of homogenous turbulence is that we can 
separate the system of Fourier modes identified by k G /C into two mutually noninteracting 
sets and K but when we do so, only the set KPl contains dynamically interacting 
modes, while the set K . ^ contains no dynamically interacting modes at all (the modes 
associated with JC ^ are frozen). By restricting the dynamically interacting modes in a finite 
model of homogeneous turbulence initially to k G KP\ we can then allow the restricted 
system to ‘expand’ into a larger system by turning on the interaction between modes in K ^ 
and KP\ so that the final set of interacting modes is k G 1C ^ U KP = K. 

We can assign an entropy to each of the isolated systems and KP\ which have Af^ 
and Af^ modes, respectively, with Af^ = A The system is actually a system of 
A A 1 ) noninteracting subsystems, each containing only one mode. Since these modes do not 
evolve, the number of states available to each isolated mode is n(k) = 1; the number of 
states available to K ^ is no more than the product Flkexn) n (k) = 1- The entropy of each 
mode is lnn(k), so that the entropy of the whole system of modes /C^ is zero (this is 
Nernst’s theorem ). However, the value of the energy (and other invariant values) locked in 
these frozen modes is not necessarily zero. 

The entropy of K P' ) can be found using Table 5.1. Recall that the entropy is the global 
minimum of the entropy functional in Table 5.1 and that the summations there are now over 
all Af^ interacting wave vectors k G /C^ . The entropy is 

5 ( 2 ) 

= min{ 8 . 7 )} 

= a(K l 2 >,a< 2 >,/?< 2 ',7 <2) ), (5.84) 


where 


a(/C,a,p, 7 ) = 


1 

2 


c 


Af(l + ln7r) — ^2 ln^(A;, a, /?, 7 ) . 

kex 


(5.85) 
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(Here, we have used A f — 2A f .) If we add the two systems together to form K = K^ U KS 2 \ 
then the total number of k 6 K is Af = A + Af^. 

Assume that the energy (and other invariant values) locked in the frozen modes JC ^ 
is proportional to those associated with the interacting modes /C 2 ); let the proportionality 
factor be r ( e.g ., E b) = rE^). The ‘temperatures’ (a^ 2 )) -1 , (,0^) _1 , and (7^) _1 are given 
in the first section of this chapter; the pertinent expressions indicate that when expansion 
from /C (2 ) to K occurs, the temperatures will drop by approximately (1 + r)/2. Therefore, 
we see from Table 5.1 that 

<p(k, a, /3, 7 ) = : ^</>(fc,ar ( 2 ) /3 ( 2 ) ,7 (2) ). (5.86) 

Expansion of K^> into K will thus lead to the following change in entropy: 

S = 2S {2) -Af (2) ln[2/(l + r)], r > 0. (5.87) 

Clearly, if r > 1, then entropy has increased. However, for 0 < r < 1, it is unclear whether 
or not entropy might decrease. That it never decreases can be shown in a more general 
argument. 

5.3.2 Proof of the Second Law 

Here, we extend the proof of [Khinchin 49] to the case of systems with more than one tem- 
perature. The value of the entropy for a given isolated subsystem Kj of modes is determined 
by finding the global minimum of the entropy functionals given in Table 5.1: 

Sj = a{K,j,OLj^j, 7j). (5.88) 

Here, aj, (3j, and 7 j are the inverse temperatures that minimize o in (5.88), for any subsystem 
j. Let us now unite the subsystems j = a and j = b to form the larger subsystem of 
interacting modes K a i , = K a U /C5. 

The previously isolated subsystems K a and K b had entropies S a and Sb- 

Sa — @{Kai Oi a , /3 a , 7a) 

S b = < j(K b ,a b J b ,%)- (5-89) 

This new subsystem K ab will have its own inverse temperatures a ab , f3 a b, and 7 a& , and entropy 

Sab- 


Sab = @(Kabi &aby Pabi 'Jab) 

= O ( Kq , Ot a b, P ab • 7 ab) T @{Kbi (%ab) Pab : ''fab) 
^ <j(AT a , a a , Pai la) T a(K b , 0 : 5 , Pbi lb)- 


(5.90) 


70 


CHAPTER 5. THERMODYNAMICS 


The second line above occurs because a is a sum over independent k and can be partitioned 
as desired, while the third line occurs because a a tPatla and a b ,P b , % minimize the entropy 
functionals related to JC a and JC b , respectively: 

&at Pat 7a) — ®abt Pabt 'Jab) 

a(K, b ,a b ,p b ,j b ) < o(K. b , a ab , f3 ab , ^ab)- (5.91) 

In looking at (5.89), it is evident that (5.90) becomes 

S a b > S a + S b . (5.92) 

This is a general result and applies to all the cases of ideal homogeneous fluid and mag- 
netofluid turbulence considered in the present work. 

Finally, one point that is often missed needs to be emphasized. The entropy of a closed 
or quasi-closed ( i.e canonical) system has a fixed value, essentially S = lnT, where T is the 
number of ‘states’ available to the system in question. The entropy of such a system does 
not ‘increase to the largest possible value if the system is left alone,’ because once the system 
is ‘left alone,’ its entropy is fixed: the entropy is a function only of the number of accessible 
states, not a function of the instantaneous location of the system point in phase space. Only 
when two previously isolated systems, each with associated entropies (i.e., measures of the 
number of available states of each isolated system), are united into a greater system, can 
total entropy actually increase, as shown in (5.92) above. 


5.4 References for Further Reading 

Two of the best references on statistical mechanics and thermodynamics are [Khinchin 49] 
and [Landau 80]. 

An early development of the concept of entropy in computer models of ideal turbulence, 
similar to that presented here, was given in [Shebalin 82], and a more recent discussion 
is given in [Shebalin 96] (where it is shown that the entropy function of [Carnevale 81, 
Carnevale 82] is not generally commensurate with absolute equilibrium ensemble theory). 


Chapter 6 

Numerical Experiments 


The basis of a statistical theory of ideal homogeneous turbulence has been given in the 
preceding chapters. How well this theory works requires a comparison ■with results drawn 
from numerical experiments. We must use numerical experiments, rather than physical 
experiments, because ideal homogeneous turbulence is not physically realizable - no physical 
fluid or magnetofluid experiment has zero dissipation throughout its experimental sample. 
However, ideal turbulence is a model system whose equations differ from real ones only in 
the absence of a linear dissipative term (there is still a linear term in ideal MHD when a 
nonzero mean magnetic field is present). Although this linear dissipative term is critical, the 
remaining nonlinear terms are identical for real and ideal sets of equations. This motivates 
our desire for a comparison: A robust statistical theory of ideal turbulence may shed some 
light on how to obtain a corresponding theory for real turbulence, which has thus far been 
elusive. 

Here we will present numerical results concerning all of the various cases of ideal homoge- 
neous turbulence: Euler and ideal MHD turbulence, for both 2-D and 3-D. These numerical 
simulations were done on various supercomputers, and the total central processing unit 
(cpu) time used to generate the results to be described here was about 560 hours. This is 
a substantial investment and calls for a clear presentation of the methods of research, the 
results obtained, and the conclusions deduced; this is a prime motivator for this publication. 
Numerical results will be presented in this chapter, and an explanation of the somewhat 
surprising results will be given in the next chapter. 


6.1 Numerical Method 

The equations of motion to be solved are (3.44) and (3.45) for 2-D, and (3.15) and (3.16) 
for 3-D. The general technique used to numerically solve these equations is called a ‘Fourier 
spectral transform method.’ (‘Spectral’ refers to the coefficients, or ‘spectrum,’ of a Fourier 
expansion.) 
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The modal equations of motion have the form 

dU< dt' t} = A&M; k )- (6-1) 

The term F(U ; k) on the right side is nonlinear (except for a linear part that occurs in the 
MHD cases for B 0 ^ 0). In what follows, when U or U(t) is used without a wave vector k 
in the argument, it signifies the set of coefficients U(k, t). 

As an example, for the 2-D Euler case, we have U( k, t ) = <D(k, t) and 


p+q=k 

i Y z 




( 6 . 2 ) 


F(a)(t);k) = i ^ z-pxq 

k, p,q£5 

Evaluating terms such as this in k-space takes an inordinate amount of time - there are 
0(N 2 ) coefficients <D(k) and each summation has 0(N 4 ) terms, so that the total floating point 
operation count is 0(N 6 ). Instead, the quantities ipip(p) and iq£>(q) are first transformed 
by a fast Fourier transform (FFT) into the x-space quantities V^(x) and Vo;(x), then the 
product v(x) = V^(x) x Vcu(x) is formed. Then we transform back by FFT to k-space: 
v(x) — » v(k) and form the result F(U] k) = ik ■ v(k). 

In 2-D, the FFTs take 0[(A/'ln A") 2 ] operations and the x-space nonlinear product takes 
0(N 2 ) operations, so that the whole transform method takes 0[(Af 2 In iV) 2 ] operations (where 
N is the number of grid points in each dimension). The comparative efficiency between work- 
ing solely in k-space and using a transform technique is thus 0(N 6 ) versus 0[(iV 2 In iV) 2 ]: 
using FFTs reduces the number of floating point operations by 0[(Af -1 In N) 2 ]. (The reduc- 
tion is 0[(iV -1 In N) 3 ] in the 3-D cases.) Thus, while a direct evaluation of the (convolution) 
sum (6.2) may be impractical for a numerical simulation, the transform method is not, and 
in fact makes such simulations possible [Orszag 72], 

The nonlinearity of (6.2) can introduce aliasing errors if care is not taken. Here, de-aliased 
nonlinear terms were produced using the Patterson-Orszag shifted-grid method [Patterson 71], 
which entails doing twice as many FFTs and setting k max = y/2N/S (this value of k max re- 
moves all but the single-aliasing errors). Specifically, for the 2-D Euler example, 


v(p,a) 

= 

ip^(p,t)e ap 

w(q, a) 

— 

iqcD(q, t)e aq 

v(p,a) 

- FFT -4 

v (x, a) 

w(p,a) 

- FFT ^ 

w(x, a) 

G(U ; k, a) 

«- FFT - 

v(x, a) x w(x, a) 


The transform cycle in (6.3) is done twice, with a = 0 and a 
de-aliased nonlinear product is then found by: 


(6.3) 

(tt/N)(1, 1, 1). The 


F(U ; k) = - * z • [G(U\ k, 0) + e~ ih k G (U; k, h)] kF(k) 


W (k) 


1, k < k 
0, k > k 


max 

max 


(6.4) 
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which removes the remaining single- aliasing errors. (The function W (k) ensures that the 
values of all coefficients with k > k max remain zero during the simulation.) 

There is another method, called the ‘2/3’ method, that requires only a single transform 
cycle by setting k max — N/3. However, the number of de-aliased modes retained for a given 
N in the Patterson-Orszag method is twice as many as in the ‘2/3’ method for 2-D and 
2 a/ 2 = 2.83 as much for 3-D simulations. Although the Patterson-Orzsag method requires 
that two FFT tramsform cycles be done, which doubles the time needed to determine the 
nonlinear terms (6.2), it makes up for this by allowing at least twice as many modes to be 
retained in the simulation for a given value of N. 

In general, the greatest portion of time in the numerical integration from one time-step 
to another goes into the evaluation of the nonlinear terms such as (6.2). In this regard, a 
time-stepping method was chosen so that the number of evaluations of (6.2) was minimized 
per time-step. For ease of coding and to minimize required core memory, lower-order time- 
integration schemes were also used. 

The two time-integration methods used to solve the modal equations (6.1) were a second- 
order Runge-Kutta method (RK2) [Potter 73] and a third-order ‘partially corrected’ Adams- 
Bashforth method [Gadzag 76]. Time-step size At was fixed, so that after n time-steps, 
t = t n = nAf; also, t n+ 1 = t n + At and t n+ 1 / 2 = t n + At/ 2 . In addition, let us use the notation 
u n = U(k,tn), fn = F(U(t n y, k), u n+1/2 = U(k, t„+ 1/2 ) , and f n+1/2 = F(U(t n+ i /2 );k). The 
RK2 method is 


Un+1/2 

Un+1 


Un + 



U n + At/n+i/2- 


(6.5) 


The RK2 method was used in the 2-D simulations and in a minority of 3-D simulations to 
be described shortly. The RK2 method is ‘self starting,’ i.e., the necessary nonlinear terms 
f n and fn+1/2 are computed at each time-step. 

The AB3 method used has two parts: The first is an Adams-Bashforth predictor: 


Un+\ — Un + 12 16/„_i T 5 jGt. — 2 J • 

The second part is an Adams-Moulton corrector: 

fn+ 1 F(u n+ 1 , k) 

Un+1 = Un + 12 [^/n+1 + 8/n — fn - 1 • 


( 6 . 6 ) 


(6.7) 


Here, the f n+ 1 are determined from the predicted coefficients u n +i, found in (6.6), in exactly 
the same way that the nonlinear terms f n + 1/2 are determined from the coefficients fi n+ 1/2 in 
RK2: by using (6.3) and (6.4). The AB3 method is started at t = 0 by assigning 

/- 2 - /-I = fo = F(u 0 -, k). 


( 6 . 8 ) 
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Thus, AB3 is effectively started with a forward Euler method. 

The time-step size was chosen for both the 2-D (N = 32) and 3-D (N = 16) runs to be 
At = 0.001. Initial checks of the 2-D and 3-D codes were performed by explicitly checking 
the values of a few selected coefficients for several time-steps, and the long-term behavior 
of the codes was confirmed by the constancy of the integral invariants (to within canonical 
fluctuations). A comparison between ensemble prediction and numerical determination of 
the means and variances of the Fourier coefficients was thought to be another good test of 
code viability, but this had the implicit assumption that ideal homogeneous turbulence was 
ergodic. As the following discussion of results will show, this assumption was not correct. 


6.2 2-D Simulations 

We begin with 2-D simulations, run either on a Control Data Corporation Cyber 205 at 
NASA Langley or on one of the Cray 2 supercomputers then at NASA Ames [Shebalin 89]. 
These simulations were performed on a 32 2 grid, with k m i n = 1 and k max = 15.08; again, 
de-aliased nonlinear terms were produced using a shifted-grid method [Patterson 71]. Since 
k max — 15.08, there are 227 independent modes for a 2-D Euler simulation, while there are 
twice this, 454 modes, for a 2-D ideal MHD simulation. Thus, the phase space for these 2-D 
Euler runs has 227 dimensions and for the 2-D ideal MHD runs has 454 dimensions. 

The modal spectra for the various runs were initialized so that Ek = Em = 0.5. The 
modes were assigned random initial phase and satisfied 

|u(k )| 2 = |b(k )| 2 ~ k A exp(—2k 2 / ko) . (6.9) 

Here, ko = 2 was used, although the exact value of ko did not matter, as the system quickly 
evolved to ‘thermodynamic equilibrium.’ 

Time-integration was performed with the RK2 method given in (6.5) and time-step size 
was At = 0.001. (The Euler runs were approximately 0.13 seconds/At, while the MHD runs 
were approximately 0.38 seconds/At.) A number of short runs were done, from t = 0 to t = 
10 (10 4 time-steps), for preliminary analysis and to allow transients to subside. After this, 
three runs were continued from t = 10 to t = 510 (5 x 10 5 time-steps): NSP (an Euler run), 
CAI (ideal MHD, B 0 = 0), and CBK (ideal MHD, B 0 = B 0 x, with B 0 — 1). 

In order to determine expectation values for the modal spectra of two of these runs, the 
entropy functionals defined in Table 5.1 were minimized with respect to a variable parameter 
(&' for run NSP and S' M for CAI). At the minimum, T' = = (T) and £' M = Em — {Em)- 

Finding L/ for run NSP allowed the inverse temperatures a and /3 to be found through (5.10) 
and (5.11), respectively; in turn, finding Em for run CAI allowed a , /?, and 7 to be found 
through (5.40), (5.41) and (5.42), respectively. In the case of run CBK, no minimization 
is required because Em — Ej 2 — (E) /2, leading to 7 = 0, as is known beforehand. In 
Figure 6.1 we show the entropy functionals cr(iF) and o(R), corresponding to NSP and CAI, 
respectively, where R = £ M /{£ — Em)- In this Figure, we have reduced the number of 
variables in the arguments of the functionals to display only the essential parameter. The 
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minima are & = 0.057397 and R = 1.0758 (£ M = 0.52178), for NSP and CAI, respectively, 
while a priori we know R = 1 for CBK. 

The ‘invariants’ for these runs are: £ and Q for NSP; £, Re, and A for CAI; £ and Rc 
for CBK. The means and standard deviations (std. dev.) with respect to time for the various 
quantities are shown in Table 6.1. 


NSP 

Mean 
Std. dev. 

CAI 

Mean 
Std. dev. 

CBK 

Mean 
Std. dev. 


£ 

0.50003 

0.00001 

£ 

1.0068 

0.0043 

£ 

1.0015 

0.0009 


Q 

29.335 

0.0020 

Rc 

0.23996 

0.00002 

Rc 

0.0040135 

0.0000045 


A 

0.50000 

0.00000 

A 

0.013076 

0.001757 


T 

0.060871 

0.012804 

R 

1.0865 

0.0677 

R 

1.0025 

0.0736 


Table 6.1: Parameters for 2-D runs, from £ = 10 to £ = 510. The predicted values are & = 
0.057397 for NSP, R = 1.0758 for CAI and R = 1 for CBK. 


It is clear from this Table that the ‘invariants’ are invariant to within small fluctuations, 
while the noninvariant quantities in each case fluctuate considerably (percent fluctuation = 
100%xstd. dev. /mean). In particular, the invariants fluctuate less than about 0.1%, while 
<p and R fluctuate about 20% and 7%, respectively. Also, A has 0% fluctuation for CAI, 
while it fluctuates about 13% for run CBK, where it is not an invariant. 

We next compare the modal kinetic and magnetic energies determined by numerical 
simulation with their canonical ensemble predictions. To accomplish this, the time-averages 
and standard deviations of all modes in the runs NSP, CAI, and CBK were determined from 
£i = 10 up to £2 = 510, for a maximum period of r = £2 — £1 = 500. In these 2-D runs, the 
primary modes were those of scalar vorticity cD(k) and magnetic scalar potential a(k). 

Over a period of time r, the averages (D(k)) T and standard deviations (<j(k)) std)T of 
the vorticity coefficients are 


H k ))<U, 9l T 

<*00>L,r 


1 ri 2 

/ w(k) dt 
T Jti 

1 F lii(k) - (<0(k)) w 

T 


dt. 


( 6 . 10 ) 

( 6 . 11 ) 


Similarly, the averages (a(k)) T and standard deviations (a(k)) strfT with respect to time 
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Figure 6.1: Entropy i 
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of the magnetic potential coefficients are 

(«( k ))a W -If «( k ) dt ( 6 - 12 ) 

T Jt\ 

<“< k ))L,T = l r |fi(k) - <“( k ))«»9,r| 2 <«• (6-13) 

T Jt\ 1 

Note also that the averages and standard deviations (squared) defined here for each coeffi- 
cient can obviously be split into averages and standard deviations (squared) of the real and 
imaginary parts of the coefficients individually, if so desired. 

The average E^ 9,T ( k), random E T ^ n ( k), and total E t ^ t ( k) kinetic energies associated with 
each mode over a given period of time r are 


ET' T i k) = ^\{u(k)) avg J 

(6.14) 

E r K an ’ T ( k ) = fc- 2 (^( k )>L,r 

(6.15) 

Ef' r { k) = Ek 9,t (k) + E™ n,T (k) . 

(6,16) 

Similarly, the average (or coherent ) E^ 9,T (k), random £^“ n (k), and total £^(k) magnetic 
energies associated with each mode over a period r are 

ET T ( k) - k 2 \(a(k)) avg J 

(6.17) 

ET’ T ( k) = fc 2 (a(k))L lT>T 

(6.18) 

E\f T ( k) = ^’ T (k) + ^r T (k)- 

(6.19) 

Adding the kinetic and magnetic parts gives the modal total energy: 


E av9 ’ T ( k) = Fr T (k) + C’ T W 

(6.20) 

E ran ’ T ( k) - E™ n ' T (k) + E r M n ’ T (k) 

(6.21) 

E tot,T (k) = Ef' T (k) + E^' T {k). 

(6.22) 

(Here we have used the relations (w(k)) sW>T = | (w(k)) atd)T | and (a(k)) afat>T = | (a(k)) stdr |.) 

We can now plot these time-averaged modal energies versus the 2-D vector k and compare 
them to each other and to canonical ensemble predictions. This is done in Figures 6.2, 6.3, 
6.4, and 6.5. Figure 6.2 consists of modal total energies (6.20), (6.21), and (6.22) averaged 


from t\ — 10 to t 2 = 60 (r = 50) for runs NSP, CAI, and CBK. Figure 6.3 consists of modal 
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total ( i.e ., kinetic) energies for run NSP averaged from t\ = 10 to t 2 = 510 (r = 500). 
Figures 6.2 and 6.3 have four parts: The first three represent numerical simulation results: 
1) E av9,T (k), 2) E ran,T (k), and 3) E tot,T ( k). Number 4 is the ensemble prediction for modal 
total energy E tot,T { k); this is purely random, since according to canonical ensemble theory 
E av9 ’ T ( k) = 0, so that E tot ’ T { k) = £ ran > T (k). 

In Figure 6.2 (r = 50), some of the modal average energies seem rather large. In looking 
at Figure 6.3 (r = 500), we see that, for run NSP at least, those modal average energies 
which appeared large at r = 50 have become considerably smaller. Also, Figure 6.2 (r = 
50) clearly shows that for run CBK ( B 0 = 1) it is only those modes with k x = 0 that have 
appreciable modal average energies. This is due to the presence of Alfven waves, which 
make the phase ip of the Fourier coefficients continually sweep through all possible values 
0 < < 27 t with the Alfven frequency k x B 0 = k x for each mode corresponding to k. (Alfven 

waves are discussed in subsection 3.2.3.) 

Figures 6.4 and 6.5 consist of modal kinetic and magnetic energies averaged from t\ = 
10 to t 2 — 510 (r = 500) for runs CAI and CBK. Figures 6.4 and 6.5 have four parts. The 
first three represent numerical simulation results: 1) E^ 9 ’ T ( k) or E < ^ 9,T ( k); 2) E™ n ' T ( k) or 
J E'^ n ’ T (k); and 3) E t £ t,T ( k) or £'^’ T (k). Number 4 is the ensemble prediction for £^“ n,T (k) 
_ ^j^’^k); again, this is purely random according to canonical ensemble theory. 

In looking at Figure 6.3 (run NSP, r = 500), we see that the time average part of the 
Fourier coefficients has become much smaller than what they were in Figure 6.2 (run NSP, r 
= 50). In Figure 6.3 (run NSP, r = 500), it is also clear that the random part produced by 
the numerical experiment is quite similar to the ensemble prediction. However, in looking at 
Figures 6.4 and 6.5, corresponding to the two ideal MHD runs CAI and CBK, respectively, 
we see that there is a significant amount of energy in E^ 9,T ( k) and E^ 9,r (k), particularly 
at the lower |k| values. In fact, these have not changed much from Figure 6.2 (runs CAI and 
CBK, r = 50). Since the canonical ensemble prediction is that E a ^ 9 ' T {k) and E^ 9,T {k) are 
identically zero, something is amiss. 

It may be asked whether or not the time allowed for the numerical simulations to run 
was long enough. To answer this question, we look at the evolution of the total coherent 
energy: 


E C k(t) = (k) 

(6.23) 

k€/C' 


E C m (t) = Y. KI’' T (k)- 

ke/c ; 

(6.24) 


For runs NSP, CAI, and CBK, the coherent energies, E^(r) and E^ 4 (t), have been averaged 
from r = 50 to r = 500, in steps of 6r — 50. If these running averages do not appear to be 
converging to some constant value, then we have not run the simulations long enough; if they 
do appear to converge to some constant value, then we may believe that we have averaged 
long enough. The running averages of the coherent energies E^(t) and E^ir), with respect 
to averaging time r, is shown in Figure 6.6. 
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12 3 4 

Figure 6.3: Run NSP - 2-D Euler time-averaged modal kinetic energies: 1 average, 2 random, 
3 total; number 4 is the ensemble prediction, which is purely random; averaging time: r = 
500. 


It is clear in Figure 6.6 that simulations NSP, CAI, and CBK have been run long enough 
so that the stationary values of the coherent energies are apparent. In the case of run NSP, 
the coherent energy E%{t) appears to be converging monotonically to zero. In this 2-D Euler 
simulation, time-averages appear to match the ensemble prediction. However, in both of the 
runs CAI and CBK, the coherent energies E^(t) and E^(t) are not converging to zero, 
but instead to some significantly nonzero values. Thus, for runs CAI and CBK, numerical 
time averages do not appear to match the ensemble prediction. In these 2-D ideal MHD 
simulations, we must therefore conclude that we have evidence for nonergodicity in ideal 
homogeneous, MHD turbulence. 

To study the implications of this more closely, we can examine the behavior of single 
Fourier coefficients. A Fourier coefficient has a real and an imaginary part, and we can plot 
one versus the other from ti = 10 to £2 = 510; this produces a ‘random walk’ on a 2-D phase 
surface. In the case of run NSP, the resulting plot is a projection of the 227-dimensional 
phase trajectory onto a 2-D plane, while for runs CAI and CBK, the resulting plot is a 
projection of the 454-dimensional phase trajectory onto a 2-D plane. Here, we will look at 
only the ideal MHD runs CAI and CBK, since it is in these runs that nonergodicity appears 
to be manifested. 

For run CAI, Figure 6.7 shows the evolution of the real versus the imaginary parts of 
o)(k) for k = (1,0) and k = (0,1), while Figure 6.8 shows the evolution the real versus the 
imaginary parts of j(k) — a(k), also for k = (1,0) and k = (0,1). Analogous plots for run 
CBK are given in Figures 6.9 and 6.10. (Absolute scales in these Figures are omitted as 
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Figure 6.4: Run CAI - 2-D ideal MHD (B 0 = 0) time-averaged energies: 1 average. 2 
random, 3 total; number 4 is the ensemble prediction, which is purely random; top: kinetic, 
bottom: magnetic; the scale is different between top and bottom; averaging time: r = 500. 
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Figure 6.5: Run CBK - 2-D ideal MHD ( B 0 = 1) time-averaged energies: 1 average, 2 
random, 3 total; number 4 is the ensemble prediction, which is purely random; top: kinetic, 
bottom: magnetic; the scale is different between top and bottom; averaging time: r = 500. 
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Figure 6.6: Coherent energies for three cases of 2-D ideal turbulence: NSP (Euler), CAI 
(MHD, B 0 = 0), and CBK (MHD, B 0 = 1) - Kin: kinetic; Mag: magnetic. 


being unimportant to the discussion at hand.) These four Figures clearly show that the 
phase trajectories for ideal MHD, with either B m = 0 or B m = 1, are not centered on the 
origin in phase space, as predicted by canonical ensemble theory. This gives further concrete 
evidence for nonergodicity in 2-D ideal MHD turbulence. 

The main difference between runs CAI and CBK is that CAI has no mean magnetic field 
present ( B 0 = 0), while CBK does have a mean field ( B 0 = 1). In the discussion in subsection 
3.2.3, it was shown that if B 0 = B 0 Sc, with B 0 ^ 0, then Alfven waves are present for those 
modes with k x ^ 0. In the limit of large B 0 , equation (3.54) describes the modal dynamics 
and the appearance of the trajectories in Figures 6.9 and 6.10 would be purely circular. For 
moderate values of B 0 , we would expect that nonlinear modal interactions would modulate 
the amplitude of the Alfven waves, and that is what we see in Figures 6.9 and 6.10, for those 
modes with B 0 / 0. These ‘nonlinear Alfven waves’ appear for all modes with k x ^ 0 when 
B 0 ^ 0, and not in those modes with k x — 0, as indicated by Figures 6.9 and 6.10. 

The presence of nonlinear Alfven waves means that the associated modal averages will 
be zero. This explains the series of spikes for run CBK along the k y axis in Part 1 of Figure 
6.2, which gives the average modal energies for r = 50. In looking at Figure 6.5, it is clear 
that these spikes are still there along the k y axis, and nowhere else. (The perspective of 
Figure 6.5 is essentially down the k y axis; the perspective in Figure 6.2 more clearly shows 
the ‘fence of spikes’ along the k y axis.) 
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CO(0,1) 




Figure 6.7: Run CAI: 2-D ideal MHD ( B 0 = 0), projection of phase trajectory: top: o)r( 0, 1) 
vs u>/( 0, 1), bottom: u R ( 1, 0) vs u)/(l, 0). 
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Figure 6.8: Run CAI: 2-D ideal MHD ( B 0 = 0), projection of phase trajectory; top: a#(0, 1) 
vs d/(0, 1), bottom: a R ( 1, 0) vs a/(l, 0). 
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6.3 3-D Simulations 

A series of 3-D simulations was also undertaken, run on a Cray YMP at NASA Lang- 
ley [Shebalin 94]. These simulations were performed on a 16 3 grid, with k min = 1 and 
k max = 7.542; de-aliased nonlinear terms were again produced using a shifted-grid method 
[Patterson 71]. Since k max = 7.542, there are 428 independent modes for a 3-D Euler simula- 
tion, while there are twice this, 856 modes, for a 3-D ideal MHD simulation. The phase space 
for 3-D Euler runs thus has 428 dimensions and for 3-D ideal MHD runs has 856 dimensions. 
The initial condition were as described by (6.9). 

Time-integration was performed with the AB3 method, given in (6.6) and (6.7), and 
time-step size was A t = 0.001; a few cases were rerun with the RK2 method (6.5), also with 
A t = 0.001; the results were essentially the same. Although AB3 required slightly more 
run-time memory than RK2, the computer execution time per At was about half as long for 
AB3 as for RK2: 0.12 sec/A t vs. 0.25 sec/A t for the Euler runs, and 0.35 sec/ A t vs. 0.71 
sec/A t for the ideal MHD runs. Overall, AB3 seems to be the more efficient method. 

Four ‘short’ runs were done, from t = 0 to t = 250 (2.5 x 10 5 time-steps), and four ‘long’ 
runs from t = 0 to at least t = 750 (7.5 x 10 5 time-steps). The runs are as follows (at the 
end of the run identifiers, ‘A’ and ‘R’ stand for AB3 and RK2, respectively). The short runs, 
which will not be discussed in much detail here, were: El A (Euler, Hk ~ 0), E2R (Euler), 
MIR (ideal MHD, B 0 = 0), and M3A (ideal MHD, B 0 — 0, k m i n — 2). The long runs were: 
E2A (Euler), MIA (ideal MHD, B 0 = 0), M2A (ideal MHD, B 0 = 0, small H M ), and M4A 
(ideal MHD, B 0 = H 0 x, with B 0 = 1). The characteristics of these runs are shown in Table 
6.2. The time-step size (At = 0.001) was sufficiently small so that the integral invariant 
values listed in Table 6.2 varied only by a few parts per million during any of the runs. 


Run 

k 2 

n 'min 

Bo 

Time 

£ 

Hr 

H c 

Hm 

E1A 

l 


250 

0.5000 

-0.04443 



E2A 

i 


750 

0.5000 

3.142 



E2R 

l 


250 

0.5000 

3.142 



MIA 

l 

0 

1000 

1.0000 


0.1326 

0.2129 

MIR 

l 

0 

250 

1.0000 


0.1326 

0.2129 

M2A 

l 

0 

750 

1.0000 


0.1326 

0.09973 

M3A 

2 

0 

250 

0.9983 


0.1258 

0.09899 

M4A 

1 

1 

750 

1.0000 

. . . 

0.1326 



Table 6.2: Parameters for the 3-D runs. 

In Table 6.2, E2A and E2R had the same initial conditions, as did runs MIA and MIR; 
these pairs of runs were done to compare the two different numerical integration methods, 
AB3 and RK2. As will be seen ahead in Tables 6.3 and 6.5, the two different methods 
produce essentially the same results (AB3 being twice as fast per time-step as RK2). Runs 
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Figure 6.11: Time evolution of magnetic helicity for 3-D ideal MHD runs M2A ( B 0 = 0) and 
M4A (B 0 = 1). 


M2A, M3A, and M4A had the same initial conditions, except that Run M3A had the k 2 =1 
coefficients set to zero; since M3 A had kN n = 2, the k 2 =1 coefficients did not evolve with 
time. Run M3 A was used to test the effect of raising k min from 1 to 2; as is evident in Tables 
6.3 and 6.5, there were some slight effects, but no major ones. 

Although runs M2A and M4A had the same initial conditions, run M4A had a mean 
magnetic field present ( B 0 = 1), while run M2A did not ( B 0 = 0). The expectation, from 
canonical ensemble theory, was that the magnetic helicity Hm would be invariant for M2A 
but not for M4A. This prediction was borne out, as is clearly seen in Figure 6.11. 

In Table 6.3 we show the inverse temperatures a, /?, and 7 associated with each of the runs 
in Table 6.2. The inverse temperatures are parameterized by 12, for 3-D Euler turbulence, 
and by £ M , or equivalently by R — £m/£k, for 3-D ideal MHD turbulence, as detailed in 
Section 5.1. Ensemble predictions and time averages were seen to be qualitatively different in 
the 2-D runs just discussed, and, as will be seen shortly, they are also qualitatively different 
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in the 3-D runs. 

In the 3-D runs, the inverse temperatures ct, (3 and 7 were determined by minimizing the 
root-mean-square (RMS) difference between ensemble and simulation modal energies: for 
the Euler runs, Q was varied to find a value Ll RM s which produce minimization, and for the 
ideal MHD runs, R was varied to find the value Rrms which provided minimization. These 
values of LIrms and Rrms are given, along with the corresponding inverse temperatures a, 
(3 and 7, for each of the runs in Table 6.3. 


Run 

Qrms 

Rrms 

a 

P 

7 

E1A 

17.08 


2.624 

0.006825 


E2A 

23.74 


15.58 

-2.062 


E2R 

23.79 


15.41 

-2.035 


MIA 


1.506 

1.774 

-0.7825 

-1.682 

MIR 


1.506 

1.773 

-0.7824 

-1.681 

M2A 


1.206 

1.547 

- 0.7142 

- 1.451 

M3A 


1.301 

1.616 

-0.7205 

-2.131 

M4A 


1.000 

1.401 

-0.7075 



Table 6.3: Inverse temperatures for the 3-D runs. 

Although the full 3-D spectra are not easily displayed as the full 2-D spectra (see Figures 
6.2, 6.3, 6.4, and 6.5) we can still get a crude measure of how well canonical predictions and 
numerical averages compare. This is done by defining ‘directionally averaged’ power spectra 
as follows: 


Ek(k) = 1 E (|u(k)| 2 ) 

^ |k|=fc 


(6.25) 

(6.26) 


The quantity n is the number of kinematically nonzero modes with |k| = A: and the averaging 
time is the total simulation time for the run in question. Figure 6.12 shows the directionally 
averaged power spectra for runs M2A and M3A. Here we see two things: First, the value of 
k m i n is not critical ( k m i n = 1 for M2A and k m i n — V 2 for M3A). Second, there is a small but 
nonnegligible difference between canonical predictions and numerical averages, which bears 
further investigation. (In fact, it was the persistence of such differences in early 2-D MHD 
simulations [Shebalin 82, Shebalin 83] which first indicated the need for a closer look.) 

At this point, we again study the behavior of single Fourier coefficients with respect to 
time. In Figure 6.13, we plot real versus imaginary parts for each of the coefficients cDj,(l, 0, 0), 
c5 z (l,0,0), £>,,(1.0,0) and £> z (l,0,0), from t\ = 1 to £2 = 1000 for run MIA. This is an MHD 
run, so the phase space has a dimension of 856. It is clear from Figure 6.13 that nonergodicity 




92 


CHAPTER 6. NUMERICAL EXPERIMENTS 


is again present, as significant mean- values (in comparison with the standard deviations) are 
apparent. Also, in Figure 6.13, the coefficients appear to have the relationships 


4(1, 0,0) = -idi, (1,0,0) 

6 y (l,0,0) “ -4(1, 0,0). (6.27) 

We also have the approximate similarity 

4 ( 1 , 0 , 0 ) ~ by ( 1 , 0 , 0 ) 

4(1, 0,0) - Ml, 0,0). (6.28) 

In fact, all coefficients with k 2 = 1 in runs MIA and M2A ( B 0 = 0) have this strong 
relationship. A quantitative measure of the strength of these apparent relationships can be 
achieved by finding the average modal correlation, defined as 


C(v, w; A; 2 ) 


1 v 1 r v(M)-w*(M) 
M 2 ) |k ^W o |v(k,t)| |w*(k,t)| ' 


(6.29) 


Here, n{k 2 ) is the number of coefficients with a set value of k 2 , and r is the averaging time. 
We draw v and w from the set b, j, u and w; the average modal correlations, as defined by 
(6.29), for the first three values of k 2 are given in Table 6.4. 


Run 

B o 

r 

k 2 

CUb;*: 2 ) 

C( u, b; k 2 ) 

C{u,u;k 2 ) 




1 

0.990 

0.896 

0.769 

MIA 

0 

1000 

2 

0.558 

0.283 

0.0497 




3 

0.439 

0.236 

0.0258 




1 

0.954 

0.776 

0.581 

M2A 

0 

750 

2 

0.535 

0.269 

0.0417 




3 

0.427 

0.231 

0.0208 




1 

0.0180 

0.239 

0.0745 

M4A 

1 

750 

2 

0.00456 

0.215 

-0.00322 




3 

-0.0154 

0.211 

-0.0195 


Table 6.4: Average modal correlations versus k 2 for three 3-D ideal MHD runs. 

It is clear in looking at Table 6.4 that the k 2 =1 modes are strongly correlated for B 0 = 
0, but not for B 0 = 1. The average modal correlations between j and b and between u and 
b for k 2 = 2, 3 are decreasing but still substantial for B 0 — 0; this is a reflection that Hm 
and He are invariants for B 0 = 0. In contrast, the correlation for k 2 = 2, 3 between u and 
uj for B 0 = 0 is negligible (a reflection that H K is not an invariant for ideal MHD). In the 
B 0 = 1 case, only u and b appear correlated, and then only moderately; for B 0 = 1, H c is 
invariant but H M and H K are not. These results suggest that u ~ V x u and b ~ V x b. 
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i.e., the fields u and b are tending toward force- free configurations, but only for B 0 — 0 and 
only for the lowest values of k 2 . We thus have a ‘low-/c 2 , quasi-force-free flow’ as a long-term 
state of an ideal magnetofluid with no mean field. (Data allowing similar conclusions for the 
Euler runs were not gathered.) 


Run 

E1A 

E c 

0.00301 

E°m 

H c k 

-0.00110 

HS 

h& 

E2A 

0.0416 


0.305 



E2R 

0.0413 


0.303 



MIA 

0.0880 

0.179 


0.0364 

0.177 

MIR 

0.0919 

0.172 


0.0365 

0.170 

M2A 

0.00562 

0.0656 


0.0158 

0.0635 

M3A 

0.00464 

0.0542 


0.0126 

0.0368 

M4A 

0.00146 

0.00148 


0.000421 

0.0000955 


Table 6.5: Coherency for 3-D runs for r = 250. 

Coherent energies for the 3-D runs can be defined in exactly the same way as for the 2-D 
runs, i.e., by (6.23) and (6.24). Furthermore, coherent helicities can be similarly defined, 
by using the time-averaged values of the Fourier coefficients in their calculation. Table 6.5 
presents the coherent energy and helicity values for an averaging time of r = 250. Runs 
MIA and MIR, in particular, appear to have a substantial amount of coherency - about 
27% in both energy and cross helicity, and 83% in coherent magnetic helicity! 

The high levels of coherent energy at r = 250, however, may decrease as the simulations 
are allowed to run longer in time. Figure 6.14 shows values of the coherent energies for 
the runs which ran to at least t = 750. In this Figure, the coherency in runs MIA and 
M2A appears to have reached a stationary value, while the coherency of runs E2A and M4A 
appears to be on the verge of reaching asymptotic steady values. In the case of run E2A, 
the behavior of the coherent energy during the time covered by Figure 6.14 is approximately 
fit by the following formula: 

E% w 0.012 + 0.033 e _t/237 . (6.30) 

Thus, the coherent energy for run E2A seems to be asymptotically approaching E% « 0.012, 
or 2.4% of the total energy of run E2A, as t — ¥ oo. 


6.4 Discussion of Numerical Results 

Combining these results with 2-D results, it is clear that ideal turbulence has, in most cases, a 
significant amount of ‘coherent energy’ (a statistically stationary nonzero value is significant 
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<0.(1, 0,0) 




*»**>•: 




H5'»' 


» . • VUSi 


by (1.0,0) 


Figure 6.13: Run MIA: 3-D ideal MHD (B 0 = 0), projection of the phase trajectory; top: 
real vs imaginary parts of w y (l,0,0) and ( 112 ( 1 , 0 , 0 ); bottom: real vs imaginary parts of 
L(1,0,0) and ^(1,0,0). 
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Averaging Time 

Figure 6.14: Coherent energies for four cases of 3-D ideal turbulence: E2A (Euler), MIA 
and M2A (MHD, B 0 = 0), and M4A (MHD, B 0 = 1); kinetic: dashed line; magnetic: solid 
line (log = log 10 ). 
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because canonical ensemble theory predicts that the value is identically zero). In other words, 
there appear to be coherent structures present in ideal homogeneous turbulence. There is an 
apparent correlation between coherent energy and the absolute value of the other invariants 
in a given run. For example, as Table 6.2 shows, runs MIA and M2A share the same values 
of E and He, but different values of Hm- Run MIA has the higher value of Hm, and also 
more coherent energy, as Figure 6.14 shows. 

Furthermore, although run M4A has the same values of E and He as runs MIA and 
M2A, it has essentially zero magnetic helicity Hm, as shown in Figure 6.11. It also has the 
lowest amount of coherent energy of the 3-D runs, as indicated by Figure 6.14. However, the 
coherent energy of run M4A can only be locked into the k x — 0 coefficients, since all other 
coefficients can be thought of as combinations of nonlinear Alfven waves, and their averages 
naturally tend to zero. (That 3-D run M4A is tending to a nonzero value of coherent energy 
is supported by Figure 6.6, in which the analogous 2-D run, CBK, clearly has reached a low, 
but nonzero stationary value of coherent energy.) 

In all of these cases, 2-D and 3-D, the anomaly is not the run that has nonzero coherent 
energy, but rather the run that has (or is clearly headed toward) a value of zero. There 
is only one case in which canonical predictions and numerical averages appear equal and it 
is the 2-D Euler case, represented here by run NSP. In looking at Figure 6.6, we see that 
the value of coherent energy for run NSP is heading monotonically to zero. This reversal 
of expectations is the central and surprising result of the numerical work detailed in this 
chapter. The question then arises: Why are the canonical ensemble predictions not the 
same as the time averages for all cases of ideal turbulence, except for 2-D Euler turbulence? 
In other words, how does nonergodicity arise in ideal homogeneous turbulence? 

The answer to this question is given in the next chapter. 


6.5 References for Further Reading 

The numerical results cited in this chapter are drawn from [Shebalin 82], [Shebalin 83], 
[Shebalin 89], and [Shebalin 94] (primarily from the latter two references). Independent 
work on ideal 3-D MHD turbulence was done by [Stribling 90], although the B 0 = 1 case 
was not discussed and time averages were not taken. 

A discussion of spectral methods is given in [Gottlieb 77], [Canuto 88], and [Boyd 01], while 
general discussions on numerical integration of differential equations can be found in many 
sources, including [Potter 73] and [Iserles 96]. 


Chapter 7 


Broken Ergodicity 


The statistical theory and computational study of ideal homogeneous turbulence has been 
laid out in the preceding chapters. The comparison of theoretical predictions with numerical 
results presented in the last chapter indicates that there is a significant mismatch between 
the two, which challenges the assumption of ergodicity implicit in the statistical theory. In 
the current chapter, we will see why this implicit assumption is generally incorrect. In what 
follows, we will demonstrate that ideal homogeneous turbulence is non-ergodic, except in the 
2 -D Euler case, and that this situation is caused by dynamically broken symmetry. Thus, 
ideal homogeneous turbulence has, inherent in its dynamics, a broken ergodicity. 


7.1 Symmetry Under P, C, T 

The classical symmetries are parity , or space reflection, P: x -» — x; charge reversal C: 
q — >• — q\ and time reversal T: t -» — t. (In addition, we can include the identity operation 
I.) The effect of P, C and T on u, w, b and j, as well as on the helicities H K = (u, u>)/ 2 , 
H c = (u, b )/2 and H M — (a, b)/ 2 , is given in Table 7.1. In Table 7.1, ’ or £ +’ signs 

indicate that the quantity in question does or does not, respectively, change sign under P, 
C, or T. 

The energies E K = (u, u)/2 and E M = (b,b)/2, the enstrophy = (u>,u;)/2, and the 
mean-squared magnetic potential A = (a, a ) /2 are all positive-definite quantities and thus 
unaffected by P, C, and T: E K , E M (and thus E = E K + E M ), fi, and A behave like scalars 
under P, C or T. Since it will be pertinent shortly, the inverse temperatures a associated 
with E in all cases of ideal turbulence are also scalars under P, C, or T, as can be seen in 
examining the relationships (5.10), (5.25) and (5.57). Similarly, examining (5.11) and (5.42) 
reveals that the inverse temperatures /3 - associated with Q in the 2 -D Euler case - and 7 - 
associated with A in the 2 -D ideal MHD case - are also scalars under the classical symmetry 
transformations. 

The helicities, however, behave differently. As Table 7.1 shows, the helicities Hi ( i = 
K, C, M) behave like pseudoscalars under P; He also transforms like a pseudoscalar under 
C. Thus the helical invariants H K , H c , and H M of ideal homogeneous turbulence are 
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I 

P 

C 

T 

u 

+ 

~ 

+ 

— 

U) 

+ 

+ 

+ 

— 

b 

+ 

+ 

— 

— 

j, a 

+ 

— 

— 

— 

H k 

+ 

— 

+ 

+ 

He 

+ 

— 

— 

+ 

H m 

+ 

— 

+ 

+ 


Table 7.1: Effect of the classical symmetry transformations P, C, and T on the signs of 
various quantities. 


pseudoscalars under at least one of the classical symmetry transformations. If we look at 
the inverse temperatures associated with the various helical invariants, we see that they, too, 
are pseudoscalars: 


3 -D Euler, eq. (5.26): 

2 - D MHD, eq. (5.41): 

3- D MHD, eq. (5.58): 

eq. (5.59): 


a Hk 

e = ~ n a 


/? = -2* C a 

0 - 2 ?* 


7 


£ — 2 Em 

'Hm 


a 


(7.1) 


The pseudoscalar nature of the inverse temperatures, on the left sides of the equations in 
(7.1), follows because, on the right sides, the inverse temperatures a (associated with total 
energy), as well as the energies and the enstrophy, are all scalars under P, C, or T, while 
the helicities are pseudoscalars. The net result is that / 3H K , (3H C} and 7 H M are all scalars. 
Thus, P, C , and T do not affect the canonical probability densities (4.30). 

An equivalent statement is that, in forming an expectation value (see Section 4.7), all of 
phase space is averaged over, i.e., regions of phase space associated with both signs of a given 
helicity are to be found in the domain of integration. This guarantees that the canonical 
ensemble theory of ideal homogeneous turbulence is invariant under the classical symmetries 
P, C, and T. 
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7.2 Broken Symmetry 

It can readily be shown that P, C, and T do not affect the form of the evolution equations 
of ideal turbulence: 

Vx(uxu+jxB) 

V x (u x B) . (7.2) 

Thus, the evolution equations (7.2) of ideal turbulence, as well as the canonical probability 
functions, are unaffected by P, C, or T. 

However, this theoretical symmetry is not duplicated in any realization of ideal homo- 
geneous turbulence containing a helical invariant. Once the system, as represented by a 
numerical model, is given an initial condition, it is also given a fixed sign of the correspond- 
ing helical invariant. The evolving system maintains the values of all associated invariants 
and their signs to within canonical fluctuations. The symmetry in the theory is not exhibited 
dynamically, and is thus called a broken symmetry. 

In phase space, this broken symmetry is manifested in the presence of effectively disjoint 
components. Although canonical probability densities give finite probabilities to the appear- 
ance of the system point anywhere in phase space, if the invariant helicity has a magnitude 
large compared to canonical fluctuations, the probability that the system point will travel 
from a component with positive helicity to one with negative helicity is effectively zero. (In 
other words, the components are not fuzzy enough to have statistically significant overlap.) 

In those cases of ideal homogeneous turbulence that possess a helical invariant, we can 
define a set-theoretic characteristic function, which indicates disjointness. Let hi — Hi/\Hi\, 
i = K,C, M, for Hi ^ 0, and hi = 0, i = K, C, M, for Hi = 0; then define 

xt = 2 < 1 + A i) 

X.“ = ’(l-W- (7.3) 

Thus, xt = 1 if the system point is on the positive helicity component and 0 if it is not. 
Similarly, xT = 1 if the system point is on the negative helicity component and 0 if it is not. 
If hi = 0, then xt — xt> an d i n a h cases we have xt E xt = !• 

For every invariant helicity Hi 0, the characteristic functions (7.3) tell us that the 
disjointness of phase space increases by a factor of two. The number D of disjoint components 
is therefore D = 2 N , where N is the number of nonzero helical invariants present. The phase 
space structures of the different cases of ideal homogeneous turbulence are given in Table 
7.2. 

Table 7.2 now helps explain the high levels of coherent energy seen in Table 6.5 and 
Figure 6.14. The highest levels are in the 3-D ideal MHD runs, such as MIA and M2A, 
whose phase spaces have four disjoint components: MIA, in particular, has a stationary 


du 

dt 

dB 

dt 
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Case 

B 0 

Invariants 

D 

2-D Euler 

— 

E, 0 

1 

3-D Euler 

— 

e,h k 

2 

2-D MHD 

0 

E, H c , A 

2 

2-D MHD 

1 

E, He 

2 

3-D MHD 

0 

E, H c , H M 

4 

3-D MHD 

1 

E,H C 

2 


Table 7.2: Phase space structure in ideal homogeneous turbulence (Hi 0, i = K,C,M ): 
D is the number of disjoint components. 

coherent energy that is about 19% of the total energy, as indicated in Figure 6.14. Run M2A 
has a coherent energy that is about 6% of the total available energy, as also seen in Figure 
6.14. In Table 6.2, we see that the major difference between the two runs is that MIA has a 
value of Hm that is more than twice that of run M2A. Although MIA and M2 A have phase 
spaces with four disjoint components, the components in MIA may be thought of as being 
‘further apart.’ 


7.3 Coherent Structure 

Classical ergodic theory ([Khinchin 49, pp. 19-38], [Sinai 94]) states that a statistical system 
is ergodic if and only if it its phase space does not consist of disjoint components (this is 
Birkhoff’s theorem). In ideal homogeneous turbulence with helical invariants, phase space 
has, as we have just seen, either two or four disjoint components. Thus, ideal homogeneous 
turbulence with helical invariants is nonergodic. The symmetry of the governing theory, as 
we have also seen, is dynamically broken by the presence of helical invariants. The concepts 
of broken symmetry and nonergodicity are combined when ideal homogeneous turbulence 
with helical invariants is said to display broken ergodicity, a term first used by [Palmer 82]. 

When a system has broken ergodicity, there is an underlying structure of disjoint com- 
ponents in phase space. If the system point is allowed to range over all components (as it 
does in an ensemble prediction), then modal ensemble averages are, by construction, auto- 
matically zero. However, in a dynamical time-average, as found by numerical simulation, the 
system point can range over only one component, and modal time-averages are not required 
to be zero. This leads to nonzero coherent energy;, combined with the otherwise chaotic 
behavior of the Fourier modes, this indicates that these modes are random variables with 
nonzero means. The structure in either k-space or x-space associated with these nonzero 
modal means may be thought of as a coherent structure. 

Therefore, broken ergodicity is equivalent to coherent structure in ideal homogeneous 
turbulence with helical invariants. 
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7.4 References for Further Reading 

A mathematical presentation of symmetry and invariance in nonlinear dynamical systems 
can be found in [Golubitsky 88]. 
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BROKEN ERGODICITY 


Chapter 8 
Epilog 


The primary purpose of this work has been to present the statistical theory of ideal homoge- 
neous turbulence. This has been done in the preceding chapters; the principal quantitative 
result found was that ideal turbulence with a helical invariant exhibits broken ergodicity. 
Thus, the assumption of ergodicity implicit in canonical ensemble theory is incorrect, and 
any predictions of this theory are only approximate. Full correctness could be restored if 
ensemble predictions were made by integrating only over those parts of phase space that had 
a definite sign to a given helical invariant. How to do this is not clear (it may not even be 
possible) because ensemble predictions are made at the level of a mode, and modal helicity 
appears unrestricted in sign, as opposed to the sum of all modal helicities, which has a fixed 
sign. 

Nevertheless, the discovery of broken ergodicity, and its corollary, coherent structure, is 
the significant result in the theory of ideal homogeneous turbulence. The question arises: Is 
broken ergodicity manifested in real, homogeneous turbulence? 


8.1 Ideal Versus Real 


In ideal turbulence, there is no physically imposed length scale, so that numerical simulations 
on a relatively small number of grid points are sufficient, as long as there are enough modes 
to ensure statistical behavior. This is equivalent to stating that the ratio p — k max jk m i n 
need only be moderate. If we take k min = 1 , then ensemble predictions of modal values can 
be renormalized so that any reference to explicit grid size disappears. 

As an example, consider (4.40), the 2-D Euler prediction for modal vorticity. If we define 
k = k /k max , a = Qt/k max and (3 = (3, the expectation value of modal vorticity can be 
rewritten as 


(Mk)l 2 ) 


k 2 

a -I- (3k 2 


(8.1) 


All other modal ensemble predictions can be renormalized in the same way. The result is 
that we always have 0 < k < 1; increasing k max only adds more points between 0 and 1, but 
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does not change the shape of the curve (8.1). 

In real homogeneous turbulence, the situation is completely different, because the dissi- 
pation coefficients v and p (viscosity and resistivity, respectively) set the size of the smallest 
length-scales. Consider Navier-Stokes turbulence (Euler turbulence with a ‘switched-on’ vis- 
cosity, v > 0). If it is stationary and dissipative, there is a balance between energy input at 
the largest scales and dissipated at the smallest scales; as is customary, let us call the energy 
being dissipated per unit density and unit time e. If there is no energy input, then the 
turbulence merely decays from some initial state, and e = \dE/dt\. As is well known, there 
is a dissipation wave number kp ~ (s/is 3 ) 1 / 4 : since \dE/dt\ = 2nd = e, as (3.75) shows, then 
k D ~ (O/zA) 1 / 4 . It was shown in [Kolmogorov 41] that, for low viscosity, at length scales in- 
termediate to the largest and the smallest scales of a flow, the energy spectra integrated over 
direction in k-space was E(k) ~ e 2 ^k~ b ^. Although many physical experiments have borne 
this law out, none of the energy spectra for ideal homogeneous turbulence is commensurate 
with it. 


8.2 Real Simulations 

In any simulation of real, homogeneous turbulence, care must be taken that k max — ko\ this 
effectively sets the size of u. In general, the smaller v is, the larger kp is, but for a realistic 
simulation ko cannot be larger than about N/2, the largest value of k along any coordinate 
axis. Therefore, in order to reduce dissipation as much as possible in a realistic simulation, 
the grid size must be made as large as practical. However, computers are finite machines and 
at any given time, there is always a limit on the number of grid points N along a coordinate 
axis. 

Currently 3-D simulations with N = 2 10 = 10 3 or 2-D simulations with N — 2 15 = 3 x 10 3 
are possible. These require long run times and in the 3-D case only achieve a moderate 
value of dissipation wave number, kp ~ 2 9 = 512. In the 2-D case, however, we have 
ko ~ 2 14 = 16,384. Matching ko with k max in any simulation is done by trying a value 
of v, running the code for a short time and seeing what f2 and e are, then adjusting v and 
rerunning, etc., for a few iterations, until v is set so that k max ~ ko- Many simulations of 
real turbulence have been done and we will not attempt to review them all here, particularly 
since they give no evidence as to ergodicity or nonergodicity. 

However, let us mention the one set of preliminary numerical results we know of that is 
pertinent to the search for broken ergodicity in real homogeneous turbulence. These results 
come from a set of runs that followed the 2-D ideal turbulence simulations presented earlier. 
The initial conditions of the 2-D ideal runs CAI and CBK, which were detailed in Chapter 
6, were used to initiate the dissipative runs CAID and CBKD, respectively. The runs CAID 
and CBKD had the same grid and time step sizes as CAI and CBK, but the viscosity u and 
resistivity r} were ‘switched on’ and the terms —uk 2 oj( k) in (3.44) and —r)k 2 a(k) in (3.45) 
were retained in the simulations. (We used v — rj = 0.01.) 

The numerical procedure was extended by treating the dissipative terms implicitly: once 
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the iterations in (6.6) and (6.7) were performed, an additional operation was done: 

(1 + nk 2 At)~ 1 u n+ i(k) —¥ u n+ i{k) 

(1 + fj,k 2 At)~ 1 u n+ i(k) -)• £t„ + i(k). (8.2) 

Here, /x is either u or rj, as appropriate. 

The dissipative simulations CAID and CBKD were run for 200 simulation times each, 
and the maximum Kolmogorov dissipation wave number k D (for a short time during the 
beginning of each run) peaked at between 26 and 46. The energy fell rapidly, as it was 
dissipated at the higher wave numbers. However, an examination of the evolution of the 
lowest wave number coefficients D(0, 1), o)(l,0), 6(0,1) and 6(1,0) with respect to time 
shows some very interesting behavior. In Figure 8.1, the evolution of <5(0, 1) and cD(l, 0) for 
run CAID is shown (compare this with Figure 6.7 for run CAI). In Figure 8.2, the evolution 
of 6(0, 1) and 6(1, 0) for run CAID is shown (compare this with Figure 6.8 for run CAI). In 
Figure 8.3, the evolution of D(0, 1) and <£>(1,0) for run CBKD is shown (compare this with 
Figure 6.9 for run CBK). In Figure 8.4, the evolution of 6(0, 1) and 6(1,0) for run CBKD is 
shown (compare this with Figure 6.10 for run CBK). 

In general, the range of the values of w(0, 1), <5(1, 0), 6(0,1) and 6(1,0) in Figures 8.1, 
8.2, 8.3 and 8.4 are about one half the range for the same coefficients in Figures 6.7, 6.8, 6.9 
and 6.10 (see [Shebalin 89] for more detail). Qualitatively, the figures indicate that, except 
for the Alfven wave behavior of <5(1,0) and 6(1,0) shown in Figures 8.3 and 8.4, the phase 
trajectories of these 2-D simulations of real homogeneous appear nonergodic. However, the 
dissipation coefficients are relatively large and the appearance of nonergodicity is suggestive, 
rather than conclusive. 


8.3 Some Future Directions 

Although 2-D Navier-Stokes turbulence contains no ideal helical invariants, 2-D MHD does, 
and it should prove very interesting if a realistic 2-D MHD simulation similar to the simu- 
lation of [Matthaeus 91], ‘Decaying, two-dimensional, Navier-Stokes turbulence at very long 
times,’ produced some interesting results in 2-D real turbulence. This work followed 2-D 
Navier-Stokes turbulence, simulated for a relatively long time - over several hundred ‘eddy 
turnover times’ - on a 512 2 grid with a maximum k D = 230. The initial flow was small-scale 
and highly chaotic, and the novel result was that only two large vortices of opposite sign 
remained at the end of the simulation. 

In the simulation of [Matthaeus 91], relaxation to a large-scale coherent structure was 
observed. In the case of 2-D MHD, the existence of an ideal helical invariant may or may not 
have an effect on the outcome - this is an open question. At the least, running on a large grid 
size (with smaller values of v and rj) would shed some more light on the apparent nonergod- 
icity in real homogeneous MHD turbulence seen in Figures 8.1, 8.2, 8.3 and 8.4 (produced 




Figure 8.1: Run CAID: 2-D real MHD ( B 0 = 0). projection of phase trajectory; top: o>r( 0, 1) 
vs cu/(0, 1), bottom: w/j(l,0) vs u)/(l,0). 
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Figure 8.2: Run CAID: 2-D real MHD ( B 0 = 0), projection of phase trajectory; top: a R (0, 1) 
vs a/(0, 1), bottom: 6^(1, 0) vs a/(l,0). 
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by relatively large values of v and rj). A 2-D real MHD simulation is also interesting for an- 
other reason: 3-D Navier-Stokes has an ideal helical invariant, while 2-D Navier-Stokes does 
not; since 2-D MHD turbulence possesses helical invariants, 3-D Navier-Stokes turbulence is 
closer, in this sense, to the 2-D real MHD case than to the 2-D Navier-Stokes case. 

Another direction for future investigation begins with the Navier-Stokes modal equations 
(3.15) with B = 0: 


= S(u, Qj\ k) - vk 2 u>(k). (8.3) 

C Lb 

In this equation, replace u(k) and u)(k) with 

v(k) = e" fc2< u(k), w(k) = e ukH w(k). (8.4) 

The result is 

dw(k) 

df = S(v, w; k). (8.5) 

Here, the nonlinear term denoted by S is the vector convolution: 

p+q= k , 

S(v, w; k) = ikx [e _2i,/p ' qt v(p) x w(q)j . (8.6) 

k, p, q 

The nonlinear term would look like the one for Euler flow, except for the presence of the 
factor exp(— 2wp-qf) within the summation (8.6). It can easily be seen that (8.5) satisfies a 
Liouville equation because dS(k)/<9v(k) = 0. However, it is an open question as to whether 
(8.5) possesses any invariants; if it did, then a canonical theory could be built for Navier- 
Stokes turbulence, just as one was built for Euler turbulence. (This would be a first step 
toward a theory of real MHD turbulence.) 

8.4 Summary 

Geometrically, an ideal invariant defines a hypersurface in phase space; the intersection of 
these hypersurfaces contains those regions in phase space where the phase point will be 
found (to within statistical fluctuations). In the presence of helical invariants, the regions 
of intersection are disjoint - each region corresponds to one sign of the helicity (as we have 
seen, canonical ensemble theory is invariant with respect to the sign of a helical invariant, 
and so integrates over both to produce an ensemble average). Again, this is the cause 
of nonergodicity: disjoint components in phase space (Birkhoff’s theorem). The ensemble 
average thus does not match the time average (as the numerical results examined in Chapter 
6 show). 

If dissipation is turned on, these hypersurfaces start collapsing; the smaller the viscosity 
and resistivity, the slower the collapse. However, the speed of the collapse is higher at larger 
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magnitudes of the wave number k. If a simulation can be done with the largest k max and kp 
possible, then presumably the collapse will be the slowest possible. If, during this collapse, 
the components in phase space are able to maintain their disjointness fairly well, then the 
phase point may be effectively stuck on the component where it started, even if it has enough 
time, dynamically speaking, to have moved a considerable distance in phase space. 

The manifestation of nonergodicity lies in significantly nonzero average values for the 
Fourier coefficients. (Significant can be defined as a large ratio between the absolute value of 
the average and the standard deviation of a Fourier mode.) Thus, rather than seeing Fourier 
mode values exhibiting a sequence of decreasing oscillations through zero, they may oscillate 
about some nonzero value which is itself slowly decreasing to zero. Hints of this have been 
seen in 2-D real MHD simulations on small grid sizes (32 2 ) [Shebalin 89]. It should prove 
interesting to run some large grid size/low dissipation numerical experiments in 2-D real 
MHD homogeneous turbulence. 


8.5 Conclusion 

The major portion of this work has been to describe the statistical mechanics and thermo- 
dynamics of ideal homogeneous turbulence. Furthermore, this was restricted to the incom- 
pressible flow of fluids and magnetofluids; i.e., Euler and ideal MHD turbulence. Numerical 
experiments exposed what the canonical theory had hidden: a broken ergodicity. Although 
the governing equations are symmetric under the classical symmetry transformations, this 
symmetry must be dynamically broken, because a system whose motion conserves helicity 
with a positive sign cannot also conserve it with a negative sign, and vice versa. Further- 
more, a few preliminary 2-D dissipative MHD runs have suggested that nonergodicity may 
also be a factor in real homogeneous turbulence. 

Ideal homogeneous turbulence serves as a model system for real turbulence, where, how- 
ever, dissipation is critical. The model is somewhat remote from reality, however, since 
ideal ensemble predictions produce a spectrum that is qualitatively different from that of 
the Kolmogorov k~ 5 N spectrum [Kolmogorov 41], a spectrum which has been verified in 
physical experiments. Nevertheless, by teaching us that some models of turbulence contain 
broken ergodicity and coherent structure, perhaps ideal homogeneous turbulence has taught 
us something that will ultimately be useful in our attempt to understand real turbulence, as 
well as other nonlinear physical systems. 


8.6 References for Further Reading 

The text [Landau 87, pp. 129-135] discusses the Kolmogorov theory (as well as other classical 
and modern results in turbulence theory). 
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